import pandas as pd
import numpy as np
from mat4py import loadmat
import netCDF4
import matplotlib
import matplotlib as mpl
import matplotlib.pyplot as plt
from numpy import mean
import netCDF4 as nc
import nc_time_axis
import xarray as xr
import matplotlib.cm as cm
import matplotlib.colors as colors
from matplotlib.colors import LinearSegmentedColormap
import scipy.io as spio
from matplotlib.pyplot import fill_betweenx
import os
from scipy import interpolate
import matplotlib.colors as colors
import string
import warnings
import seaborn as sns
from matplotlib.legend import Legend
#creating data arrays for the observational data
dm= spio.loadmat('/Users/nuhahameed/Documents/Data for MSci/ITPAJX_Robert .mat')
mo=dm['AJX'][0,0]['mo'][0]
nl=np.zeros(12)
TA=np.zeros([12,500])
SA=np.zeros([12,500])
for i,m in enumerate(mo):
TA[m-1][np.where(~np.isnan(dm['AJX'][0,0]['T'][i,:]))]+=dm['AJX'][0,0]['T'][i,:][np.where(~np.isnan(dm['AJX'][0,0]['T'][i,:]))]
SA[m-1][np.where(~np.isnan(dm['AJX'][0,0]['S'][i,:]))]+=dm['AJX'][0,0]['S'][i,:][np.where(~np.isnan(dm['AJX'][0,0]['S'][i,:]))]
nl[m-1]+=1
TA/=nl[:,np.newaxis]
SA/=nl[:,np.newaxis]
mo=dm['ITP'][0,0]['mo'][0]
nl=np.zeros(12)
TI=np.zeros([12,500])
SI=np.zeros([12,500])
for i,m in enumerate(mo):
TI[m-1][np.where(~np.isnan(dm['ITP'][0,0]['T'][i,:]))]+=dm['ITP'][0,0]['T'][i,:][np.where(~np.isnan(dm['ITP'][0,0]['T'][i,:]))]
SI[m-1][np.where(~np.isnan(dm['ITP'][0,0]['S'][i,:]))]+=dm['ITP'][0,0]['S'][i,:][np.where(~np.isnan(dm['ITP'][0,0]['S'][i,:]))]
# TI[m-1]+=dm['ITP'][0,0]['T'][i,:]
# SI[m-1]+=dm['ITP'][0,0]['S'][i,:]
nl[m-1]+=1
TI/=nl[:,np.newaxis]
SI/=nl[:,np.newaxis]
depth_res = 1 # depth resolution in meters
depth_values = np.arange(500) * depth_res # create array of depth values
AJX_avg_salinity = np.mean(SA, axis=0) # compute average salinity across all months
ITP_avg_salinity = np.nanmean(SI, axis=0) #compute average salinity across all months
depth_range = (depth_values >= 0) & (depth_values <= 150)
#AIDJEX
AJX_salinity_depth_range = AJX_avg_salinity[depth_range] #salinity profile up to 150 meters
AJX_Phi = np.sum(AJX_salinity_depth_range - 33)
#ITP
ITP_salinity_depth_range = ITP_avg_salinity[depth_range] #salinity profile up to 150 meters
ITP_Phi = np.sum(ITP_salinity_depth_range - 33)
print("Phi for AIDJEX:", AJX_Phi)
print("Phi for ITP:", ITP_Phi)
# Calculate the ratio of Phi values for AIDJEX and ITP
Phi_ratio = ITP_Phi / AJX_Phi
print("Phi ratio (AIDJEX to ITP):", Phi_ratio)
Phi for AIDJEX: -160.5791112075099 Phi for ITP: -327.07891460261646 Phi ratio (AIDJEX to ITP): 2.0368708740699506
import matplotlib.pyplot as plt
# Extract salinity at depth 0 for each month from May to December
SA_depth0 = SA[4:12, 0] # May is at index 4 (0-based), depth 0 is at index 0
SI_depth0 = SI[4:12, 0]
# Create an array for the months (May to December)
months = ['May', 'Jun', 'Jul', 'Aug', 'Sep', 'Oct', 'Nov', 'Dec']
# Plot the data
plt.plot(months, SA_depth0, marker='o', label='AJX (SA)')
plt.plot(months, SI_depth0, marker='o', label='ITP (SI)')
# Customize the plot
plt.xlabel('Months')
plt.ylabel('Salinity at Depth 0')
plt.title('Salinity at Depth 0 (May to December)')
plt.legend()
# Show the plot
plt.show()
ratio = ITP_Phi / AJX_Phi
depth_range = (depth_values >= 0) & (depth_values <= 150)
depth_range_values = depth_values[depth_range]
# rescaled AJX profile
rescaled_AJX_salinity_depth_range = (AJX_salinity_depth_range - 33) * ratio + 33
# Plot
plt.plot(AJX_salinity_depth_range, depth_range_values, label='Original AIDJEX', color='red', linewidth=4)
plt.plot(rescaled_AJX_salinity_depth_range, depth_range_values, label='Rescaled AIDJEX', color='blue', linewidth=4)
plt.plot(ITP_salinity_depth_range, depth_range_values, label='ITP', color='orange', linewidth=4)
plt.gca().invert_yaxis()
plt.xlabel('Salinity (g/kg)', fontsize =14)
plt.ylabel('Depth (m)', fontsize =14)
plt.title('Rescaled AIDJEX Salinity Profile', fontsize=15)
plt.legend()
plt.savefig('rescaled AJX.png')
plt.show()
# Calculate derivatives for both datasets
dSI_dz_AJX = np.gradient(AJX_avg_salinity, depth_values)
dSI_dz_ITP = np.gradient(ITP_avg_salinity, depth_values)
# Find the maximum gradient and its corresponding depth for both datasets
max_gradient_AJX = np.max(dSI_dz_AJX)
max_gradient_depth_AJX = depth_values[np.argmax(dSI_dz_AJX)]
max_gradient_ITP = np.max(dSI_dz_ITP)
max_gradient_depth_ITP = depth_values[np.argmax(dSI_dz_ITP)]
# Create a dictionary to store the dataset names, maximum gradients, and corresponding depths
data = {'Dataset': ['AIDJEX', 'ITP'],
'Max Gradient': [max_gradient_AJX, max_gradient_ITP],
'Depth of Max Gradient': [max_gradient_depth_AJX, max_gradient_depth_ITP]}
# Convert the dictionary to a pandas DataFrame
df_observations = pd.DataFrame(data)
df_observations.to_excel('Observations.xlsx', index=False)
# Display the DataFrame
print(df_observations)
Dataset Max Gradient Depth of Max Gradient 0 AIDJEX 0.041562 60 1 ITP 0.126548 34
# Calculate derivatives for both datasets
dSI_dz_AJX = np.gradient(AJX_avg_salinity, depth_values)
dSI_dz_ITP = np.gradient(ITP_avg_salinity, depth_values)
dSI_rescaled_AJX_salinity_depth_range = np.gradient(rescaled_AJX_salinity_depth_range, depth_range_values)
# Plot the derivatives
plt.plot(dSI_dz_AJX, depth_values, label='1975 AIDJEX', color='red', linewidth=4)
plt.plot(dSI_dz_ITP, depth_values, label='2006-2012 ITP', color='orange', linewidth=4)
plt.plot(dSI_rescaled_AJX_salinity_depth_range, depth_range_values, label='Rescaled AIDJEX', color='blue', linewidth=4)
# Set y-axis limits and invert the y-axis
plt.ylim(150, 0)
# Add labels, a title, and a legend to the plot
plt.xlabel('Salinity Gradient (g/kg/m)', fontsize=14)
plt.ylabel('Depth (m)', fontsize=14)
plt.title('Vertical Salinity Gradient', fontsize=15)
plt.legend()
plt.savefig('vertical salinity gradiet observations')
plt.show()
import numpy as np
import matplotlib.pyplot as plt
# Assuming the 'SA' and 'SI' arrays have already been defined and filled with data
depth_res = 1 # depth resolution in meters
depth_values = np.arange(500) * depth_res # create array of depth values
max_gradient_depths_AJX = []
max_gradient_depths_ITP = []
for month_idx in range(4, 12): # Loop over months from May (index 4) to December (index 11)
# Calculate gradients for the salinity profiles of both datasets
dSA_dz_AJX = np.gradient(SA[month_idx, :], depth_values)
dSI_dz_ITP = np.gradient(SI[month_idx, :], depth_values)
# Find the depth of the maximum gradient for both datasets
max_gradient_depth_AJX = depth_values[np.argmax(dSA_dz_AJX)]
max_gradient_depth_ITP = depth_values[np.argmax(dSI_dz_ITP)]
max_gradient_depths_AJX.append(max_gradient_depth_AJX)
max_gradient_depths_ITP.append(max_gradient_depth_ITP)
# Create an array for the months (May to December)
months = ['May', 'Jun', 'Jul', 'Aug', 'Sep', 'Oct', 'Nov', 'Dec']
# Plot the depth of the maximum gradient for both datasets
plt.plot(months, max_gradient_depths_AJX, label='AJX (SA)')
plt.plot(months, max_gradient_depths_ITP, label='ITP (SI)')
# Customize the plot
plt.xlabel('Months')
plt.ylabel('Depth of Maximum Gradient')
plt.title('Depth of Maximum Gradient for Salinity Profiles (May to December)')
plt.legend()
# Show the plot
plt.show()
warnings.filterwarnings('ignore', category=np.VisibleDeprecationWarning)
input_folder ='/Volumes/Seagate Hub/cmip6_downloader/data/Salinity models/'
myfiles = os.listdir(input_folder)
for f in myfiles:
if '._so' in f:
myfiles.remove(f)
fig, axs = plt.subplots(8, 4, layout='tight', figsize=(5 * 5, 5 * 8), sharex=True)
axs.reshape(-4)
axs[-1, -1].remove()
axs[-1, -2].remove()
def int_to_roman(integer):
val = [
1000, 900, 500, 400,
100, 90, 50, 40,
10, 9, 5, 4,
1
]
syb = [
"M", "CM", "D", "CD",
"C", "XC", "L", "XL",
"X", "IX", "V", "IV",
"I"
]
roman_num = ''
i = 0
while integer > 0:
for _ in range(integer // val[i]):
roman_num += syb[i]
integer -= val[i]
i += 1
return roman_num
def figureone(file, i):
ds = xr.open_dataset(input_folder+file) #opens the dataset
if 'CESM' in file:
ds['lev'] = ds['lev']/100
model_name = file.split('so.')[-1].split('.')[0]
Year = ['1970', '1980', '1990', '2000', '2010']
for j, year in enumerate(Year):
ds_mean = ds.so.isel(time=np.arange((int(year)-1850)*12, (int(year)+1-1850)*12,1)).mean(dim='time')
y = ds_mean.mean(dim=['longitude', 'latitude']) #works out the average within the longitudanal range
ds_std = ds.so.isel(time=np.arange((int(year)-1850)*12, (int(year)+1-1850)*12,1)).std(dim='time')
ds_std_latlon = ds_std.std(dim = ['longitude', 'latitude'])
y_upper = y + ds_std_latlon
y_lower = y - ds_std_latlon
axs.reshape(-1)[i].plot(y, ds['lev'], label=year, color=plt.cm.viridis(j/len(Year)), linewidth = 4)
for year in range(1970, 2011):
ds_year = ds.so.isel(time=np.arange((year-1850)*12, (year+1-1850)*12, 1)).mean(dim='time')
y_year = ds_year.mean(dim=['longitude', 'latitude'])
axs.reshape(-1)[i].plot(y_year, ds['lev'], color='gray', alpha=0.3, linewidth=0.3)
axs.reshape(-1)[i].fill_betweenx(ds['lev'], y_upper, y_lower, color=plt.cm.viridis(j/len(Year)), alpha=0.2)
axs.reshape(-1)[i].set_xlabel('Salinity (g/kg)', fontsize = 20)
axs.reshape(-1)[i].set_ylabel('Depth (m)', fontsize = 20)
axs.reshape(-1)[i].set_ylim(0,150)
axs.reshape(-1)[i].set_xlim(27,35)
axs.reshape(-1)[i].tick_params(axis='x', labelsize=18)
axs.reshape(-1)[i].tick_params(axis='y', labelsize=18)
axs.reshape(-1)[i].xaxis.set_tick_params(which='both', labelbottom=True)
axs.reshape(-1)[i].invert_yaxis()
axs.reshape(-1)[i].set_title(model_name, fontsize = 25)
axs.reshape(-1)[i].plot(AJX_avg_salinity, depth_values, label='1975 AIDJEX', color='red', linewidth = 4)
axs.reshape(-1)[i].plot(ITP_avg_salinity, depth_values, label='2006-2012 ITP', color='orange', linewidth =4)
if i < 26:
subplot_label = string.ascii_lowercase[i]
else:
subplot_label = int_to_roman((i - 25))
axs.reshape(-1)[i].text(0.02, 1.1, subplot_label, transform=axs.reshape(-1)[i].transAxes, fontsize=16, fontweight='bold', va='top', ha='left')
i=0
for file in myfiles:
figureone(file, i)
i=i+1
Year = ['1970', '1980', '1990', '2000', '2010']
year_colors = [plt.cm.viridis(j/len(Year)) for j in range(len(Year))]
year_handles = [plt.Line2D([0], [0], color=color, lw=4, label=year) for year, color in zip(Year, year_colors)]
aidjex_handle = plt.Line2D([0], [0], color='red', lw=4, label='1975 AIDJEX')
itp_handle = plt.Line2D([0], [0], color='orange', lw=4, label='2006-2012 ITP')
legend_handles = year_handles + [aidjex_handle, itp_handle]
# Customize the position of the legend by changing x and y
x, y = 0.713, 0.12
fig.legend(legend_handles, [handle.get_label() for handle in legend_handles], fontsize=30, bbox_to_anchor=(x, y), ncol=1, borderaxespad=0)
plt.savefig('Salinity_profiles_09_04.png')
/Users/nuhahameed/opt/miniconda3/lib/python3.9/site-packages/numpy/lib/nanfunctions.py:1670: RuntimeWarning: Degrees of freedom <= 0 for slice. var = nanvar(a, axis=axis, dtype=dtype, out=out, ddof=ddof, /Users/nuhahameed/opt/miniconda3/lib/python3.9/site-packages/numpy/lib/nanfunctions.py:1670: RuntimeWarning: Degrees of freedom <= 0 for slice. var = nanvar(a, axis=axis, dtype=dtype, out=out, ddof=ddof, /Users/nuhahameed/opt/miniconda3/lib/python3.9/site-packages/numpy/lib/nanfunctions.py:1670: RuntimeWarning: Degrees of freedom <= 0 for slice. var = nanvar(a, axis=axis, dtype=dtype, out=out, ddof=ddof, /Users/nuhahameed/opt/miniconda3/lib/python3.9/site-packages/numpy/lib/nanfunctions.py:1670: RuntimeWarning: Degrees of freedom <= 0 for slice. var = nanvar(a, axis=axis, dtype=dtype, out=out, ddof=ddof, /Users/nuhahameed/opt/miniconda3/lib/python3.9/site-packages/numpy/lib/nanfunctions.py:1670: RuntimeWarning: Degrees of freedom <= 0 for slice. var = nanvar(a, axis=axis, dtype=dtype, out=out, ddof=ddof, /Users/nuhahameed/opt/miniconda3/lib/python3.9/site-packages/numpy/lib/nanfunctions.py:1670: RuntimeWarning: Degrees of freedom <= 0 for slice. var = nanvar(a, axis=axis, dtype=dtype, out=out, ddof=ddof, /Users/nuhahameed/opt/miniconda3/lib/python3.9/site-packages/numpy/lib/nanfunctions.py:1670: RuntimeWarning: Degrees of freedom <= 0 for slice. var = nanvar(a, axis=axis, dtype=dtype, out=out, ddof=ddof, /Users/nuhahameed/opt/miniconda3/lib/python3.9/site-packages/numpy/lib/nanfunctions.py:1670: RuntimeWarning: Degrees of freedom <= 0 for slice. var = nanvar(a, axis=axis, dtype=dtype, out=out, ddof=ddof, /Users/nuhahameed/opt/miniconda3/lib/python3.9/site-packages/numpy/lib/nanfunctions.py:1670: RuntimeWarning: Degrees of freedom <= 0 for slice. var = nanvar(a, axis=axis, dtype=dtype, out=out, ddof=ddof, /Users/nuhahameed/opt/miniconda3/lib/python3.9/site-packages/numpy/lib/nanfunctions.py:1670: RuntimeWarning: Degrees of freedom <= 0 for slice. var = nanvar(a, axis=axis, dtype=dtype, out=out, ddof=ddof, /Users/nuhahameed/opt/miniconda3/lib/python3.9/site-packages/numpy/lib/nanfunctions.py:1670: RuntimeWarning: Degrees of freedom <= 0 for slice. var = nanvar(a, axis=axis, dtype=dtype, out=out, ddof=ddof, /Users/nuhahameed/opt/miniconda3/lib/python3.9/site-packages/numpy/lib/nanfunctions.py:1670: RuntimeWarning: Degrees of freedom <= 0 for slice. var = nanvar(a, axis=axis, dtype=dtype, out=out, ddof=ddof, /Users/nuhahameed/opt/miniconda3/lib/python3.9/site-packages/numpy/lib/nanfunctions.py:1670: RuntimeWarning: Degrees of freedom <= 0 for slice. var = nanvar(a, axis=axis, dtype=dtype, out=out, ddof=ddof, /Users/nuhahameed/opt/miniconda3/lib/python3.9/site-packages/numpy/lib/nanfunctions.py:1670: RuntimeWarning: Degrees of freedom <= 0 for slice. var = nanvar(a, axis=axis, dtype=dtype, out=out, ddof=ddof, /Users/nuhahameed/opt/miniconda3/lib/python3.9/site-packages/numpy/lib/nanfunctions.py:1670: RuntimeWarning: Degrees of freedom <= 0 for slice. var = nanvar(a, axis=axis, dtype=dtype, out=out, ddof=ddof, /Users/nuhahameed/opt/miniconda3/lib/python3.9/site-packages/numpy/lib/nanfunctions.py:1670: RuntimeWarning: Degrees of freedom <= 0 for slice. var = nanvar(a, axis=axis, dtype=dtype, out=out, ddof=ddof, /Users/nuhahameed/opt/miniconda3/lib/python3.9/site-packages/numpy/lib/nanfunctions.py:1670: RuntimeWarning: Degrees of freedom <= 0 for slice. var = nanvar(a, axis=axis, dtype=dtype, out=out, ddof=ddof, /Users/nuhahameed/opt/miniconda3/lib/python3.9/site-packages/numpy/lib/nanfunctions.py:1670: RuntimeWarning: Degrees of freedom <= 0 for slice. var = nanvar(a, axis=axis, dtype=dtype, out=out, ddof=ddof, /Users/nuhahameed/opt/miniconda3/lib/python3.9/site-packages/numpy/lib/nanfunctions.py:1670: RuntimeWarning: Degrees of freedom <= 0 for slice. var = nanvar(a, axis=axis, dtype=dtype, out=out, ddof=ddof, /Users/nuhahameed/opt/miniconda3/lib/python3.9/site-packages/numpy/lib/nanfunctions.py:1670: RuntimeWarning: Degrees of freedom <= 0 for slice. var = nanvar(a, axis=axis, dtype=dtype, out=out, ddof=ddof, /Users/nuhahameed/opt/miniconda3/lib/python3.9/site-packages/numpy/lib/nanfunctions.py:1670: RuntimeWarning: Degrees of freedom <= 0 for slice. var = nanvar(a, axis=axis, dtype=dtype, out=out, ddof=ddof, /Users/nuhahameed/opt/miniconda3/lib/python3.9/site-packages/numpy/lib/nanfunctions.py:1670: RuntimeWarning: Degrees of freedom <= 0 for slice. var = nanvar(a, axis=axis, dtype=dtype, out=out, ddof=ddof, /Users/nuhahameed/opt/miniconda3/lib/python3.9/site-packages/numpy/lib/nanfunctions.py:1670: RuntimeWarning: Degrees of freedom <= 0 for slice. var = nanvar(a, axis=axis, dtype=dtype, out=out, ddof=ddof, /Users/nuhahameed/opt/miniconda3/lib/python3.9/site-packages/numpy/lib/nanfunctions.py:1670: RuntimeWarning: Degrees of freedom <= 0 for slice. var = nanvar(a, axis=axis, dtype=dtype, out=out, ddof=ddof, /Users/nuhahameed/opt/miniconda3/lib/python3.9/site-packages/numpy/lib/nanfunctions.py:1670: RuntimeWarning: Degrees of freedom <= 0 for slice. var = nanvar(a, axis=axis, dtype=dtype, out=out, ddof=ddof, /Users/nuhahameed/opt/miniconda3/lib/python3.9/site-packages/numpy/lib/nanfunctions.py:1670: RuntimeWarning: Degrees of freedom <= 0 for slice. var = nanvar(a, axis=axis, dtype=dtype, out=out, ddof=ddof, /Users/nuhahameed/opt/miniconda3/lib/python3.9/site-packages/numpy/lib/nanfunctions.py:1670: RuntimeWarning: Degrees of freedom <= 0 for slice. var = nanvar(a, axis=axis, dtype=dtype, out=out, ddof=ddof, /Users/nuhahameed/opt/miniconda3/lib/python3.9/site-packages/numpy/lib/nanfunctions.py:1670: RuntimeWarning: Degrees of freedom <= 0 for slice. var = nanvar(a, axis=axis, dtype=dtype, out=out, ddof=ddof, /Users/nuhahameed/opt/miniconda3/lib/python3.9/site-packages/numpy/lib/nanfunctions.py:1670: RuntimeWarning: Degrees of freedom <= 0 for slice. var = nanvar(a, axis=axis, dtype=dtype, out=out, ddof=ddof, /Users/nuhahameed/opt/miniconda3/lib/python3.9/site-packages/numpy/lib/nanfunctions.py:1670: RuntimeWarning: Degrees of freedom <= 0 for slice. var = nanvar(a, axis=axis, dtype=dtype, out=out, ddof=ddof, /Users/nuhahameed/opt/miniconda3/lib/python3.9/site-packages/numpy/lib/nanfunctions.py:1670: RuntimeWarning: Degrees of freedom <= 0 for slice. var = nanvar(a, axis=axis, dtype=dtype, out=out, ddof=ddof, /Users/nuhahameed/opt/miniconda3/lib/python3.9/site-packages/numpy/lib/nanfunctions.py:1670: RuntimeWarning: Degrees of freedom <= 0 for slice. var = nanvar(a, axis=axis, dtype=dtype, out=out, ddof=ddof, /Users/nuhahameed/opt/miniconda3/lib/python3.9/site-packages/numpy/lib/nanfunctions.py:1670: RuntimeWarning: Degrees of freedom <= 0 for slice. var = nanvar(a, axis=axis, dtype=dtype, out=out, ddof=ddof, /Users/nuhahameed/opt/miniconda3/lib/python3.9/site-packages/numpy/lib/nanfunctions.py:1670: RuntimeWarning: Degrees of freedom <= 0 for slice. var = nanvar(a, axis=axis, dtype=dtype, out=out, ddof=ddof, /Users/nuhahameed/opt/miniconda3/lib/python3.9/site-packages/numpy/lib/nanfunctions.py:1670: RuntimeWarning: Degrees of freedom <= 0 for slice. var = nanvar(a, axis=axis, dtype=dtype, out=out, ddof=ddof, /Users/nuhahameed/opt/miniconda3/lib/python3.9/site-packages/numpy/lib/nanfunctions.py:1670: RuntimeWarning: Degrees of freedom <= 0 for slice. var = nanvar(a, axis=axis, dtype=dtype, out=out, ddof=ddof, /Users/nuhahameed/opt/miniconda3/lib/python3.9/site-packages/numpy/lib/nanfunctions.py:1670: RuntimeWarning: Degrees of freedom <= 0 for slice. var = nanvar(a, axis=axis, dtype=dtype, out=out, ddof=ddof, /Users/nuhahameed/opt/miniconda3/lib/python3.9/site-packages/numpy/lib/nanfunctions.py:1670: RuntimeWarning: Degrees of freedom <= 0 for slice. var = nanvar(a, axis=axis, dtype=dtype, out=out, ddof=ddof, /Users/nuhahameed/opt/miniconda3/lib/python3.9/site-packages/numpy/lib/nanfunctions.py:1670: RuntimeWarning: Degrees of freedom <= 0 for slice. var = nanvar(a, axis=axis, dtype=dtype, out=out, ddof=ddof, /Users/nuhahameed/opt/miniconda3/lib/python3.9/site-packages/numpy/lib/nanfunctions.py:1670: RuntimeWarning: Degrees of freedom <= 0 for slice. var = nanvar(a, axis=axis, dtype=dtype, out=out, ddof=ddof, /Users/nuhahameed/opt/miniconda3/lib/python3.9/site-packages/numpy/lib/nanfunctions.py:1670: RuntimeWarning: Degrees of freedom <= 0 for slice. var = nanvar(a, axis=axis, dtype=dtype, out=out, ddof=ddof, /Users/nuhahameed/opt/miniconda3/lib/python3.9/site-packages/numpy/lib/nanfunctions.py:1670: RuntimeWarning: Degrees of freedom <= 0 for slice. var = nanvar(a, axis=axis, dtype=dtype, out=out, ddof=ddof, /Users/nuhahameed/opt/miniconda3/lib/python3.9/site-packages/numpy/lib/nanfunctions.py:1670: RuntimeWarning: Degrees of freedom <= 0 for slice. var = nanvar(a, axis=axis, dtype=dtype, out=out, ddof=ddof, /Users/nuhahameed/opt/miniconda3/lib/python3.9/site-packages/numpy/lib/nanfunctions.py:1670: RuntimeWarning: Degrees of freedom <= 0 for slice. var = nanvar(a, axis=axis, dtype=dtype, out=out, ddof=ddof, /Users/nuhahameed/opt/miniconda3/lib/python3.9/site-packages/numpy/lib/nanfunctions.py:1670: RuntimeWarning: Degrees of freedom <= 0 for slice. var = nanvar(a, axis=axis, dtype=dtype, out=out, ddof=ddof, /Users/nuhahameed/opt/miniconda3/lib/python3.9/site-packages/numpy/lib/nanfunctions.py:1670: RuntimeWarning: Degrees of freedom <= 0 for slice. var = nanvar(a, axis=axis, dtype=dtype, out=out, ddof=ddof, /Users/nuhahameed/opt/miniconda3/lib/python3.9/site-packages/numpy/lib/nanfunctions.py:1670: RuntimeWarning: Degrees of freedom <= 0 for slice. var = nanvar(a, axis=axis, dtype=dtype, out=out, ddof=ddof, /Users/nuhahameed/opt/miniconda3/lib/python3.9/site-packages/numpy/lib/nanfunctions.py:1670: RuntimeWarning: Degrees of freedom <= 0 for slice. var = nanvar(a, axis=axis, dtype=dtype, out=out, ddof=ddof, /Users/nuhahameed/opt/miniconda3/lib/python3.9/site-packages/numpy/lib/nanfunctions.py:1670: RuntimeWarning: Degrees of freedom <= 0 for slice. var = nanvar(a, axis=axis, dtype=dtype, out=out, ddof=ddof, /Users/nuhahameed/opt/miniconda3/lib/python3.9/site-packages/numpy/lib/nanfunctions.py:1670: RuntimeWarning: Degrees of freedom <= 0 for slice. var = nanvar(a, axis=axis, dtype=dtype, out=out, ddof=ddof, /Users/nuhahameed/opt/miniconda3/lib/python3.9/site-packages/numpy/lib/nanfunctions.py:1670: RuntimeWarning: Degrees of freedom <= 0 for slice. var = nanvar(a, axis=axis, dtype=dtype, out=out, ddof=ddof, /Users/nuhahameed/opt/miniconda3/lib/python3.9/site-packages/numpy/lib/nanfunctions.py:1670: RuntimeWarning: Degrees of freedom <= 0 for slice. var = nanvar(a, axis=axis, dtype=dtype, out=out, ddof=ddof, /Users/nuhahameed/opt/miniconda3/lib/python3.9/site-packages/numpy/lib/nanfunctions.py:1670: RuntimeWarning: Degrees of freedom <= 0 for slice. var = nanvar(a, axis=axis, dtype=dtype, out=out, ddof=ddof, /Users/nuhahameed/opt/miniconda3/lib/python3.9/site-packages/numpy/lib/nanfunctions.py:1670: RuntimeWarning: Degrees of freedom <= 0 for slice. var = nanvar(a, axis=axis, dtype=dtype, out=out, ddof=ddof, /Users/nuhahameed/opt/miniconda3/lib/python3.9/site-packages/numpy/lib/nanfunctions.py:1670: RuntimeWarning: Degrees of freedom <= 0 for slice. var = nanvar(a, axis=axis, dtype=dtype, out=out, ddof=ddof, /Users/nuhahameed/opt/miniconda3/lib/python3.9/site-packages/numpy/lib/nanfunctions.py:1670: RuntimeWarning: Degrees of freedom <= 0 for slice. var = nanvar(a, axis=axis, dtype=dtype, out=out, ddof=ddof, /Users/nuhahameed/opt/miniconda3/lib/python3.9/site-packages/numpy/lib/nanfunctions.py:1670: RuntimeWarning: Degrees of freedom <= 0 for slice. var = nanvar(a, axis=axis, dtype=dtype, out=out, ddof=ddof, /Users/nuhahameed/opt/miniconda3/lib/python3.9/site-packages/numpy/lib/nanfunctions.py:1670: RuntimeWarning: Degrees of freedom <= 0 for slice. var = nanvar(a, axis=axis, dtype=dtype, out=out, ddof=ddof, /Users/nuhahameed/opt/miniconda3/lib/python3.9/site-packages/numpy/lib/nanfunctions.py:1670: RuntimeWarning: Degrees of freedom <= 0 for slice. var = nanvar(a, axis=axis, dtype=dtype, out=out, ddof=ddof, /Users/nuhahameed/opt/miniconda3/lib/python3.9/site-packages/numpy/lib/nanfunctions.py:1670: RuntimeWarning: Degrees of freedom <= 0 for slice. var = nanvar(a, axis=axis, dtype=dtype, out=out, ddof=ddof, /Users/nuhahameed/opt/miniconda3/lib/python3.9/site-packages/numpy/lib/nanfunctions.py:1670: RuntimeWarning: Degrees of freedom <= 0 for slice. var = nanvar(a, axis=axis, dtype=dtype, out=out, ddof=ddof, /Users/nuhahameed/opt/miniconda3/lib/python3.9/site-packages/numpy/lib/nanfunctions.py:1670: RuntimeWarning: Degrees of freedom <= 0 for slice. var = nanvar(a, axis=axis, dtype=dtype, out=out, ddof=ddof, /Users/nuhahameed/opt/miniconda3/lib/python3.9/site-packages/numpy/lib/nanfunctions.py:1670: RuntimeWarning: Degrees of freedom <= 0 for slice. var = nanvar(a, axis=axis, dtype=dtype, out=out, ddof=ddof, /Users/nuhahameed/opt/miniconda3/lib/python3.9/site-packages/numpy/lib/nanfunctions.py:1670: RuntimeWarning: Degrees of freedom <= 0 for slice. var = nanvar(a, axis=axis, dtype=dtype, out=out, ddof=ddof, /Users/nuhahameed/opt/miniconda3/lib/python3.9/site-packages/numpy/lib/nanfunctions.py:1670: RuntimeWarning: Degrees of freedom <= 0 for slice. var = nanvar(a, axis=axis, dtype=dtype, out=out, ddof=ddof, /Users/nuhahameed/opt/miniconda3/lib/python3.9/site-packages/numpy/lib/nanfunctions.py:1670: RuntimeWarning: Degrees of freedom <= 0 for slice. var = nanvar(a, axis=axis, dtype=dtype, out=out, ddof=ddof, /Users/nuhahameed/opt/miniconda3/lib/python3.9/site-packages/numpy/lib/nanfunctions.py:1670: RuntimeWarning: Degrees of freedom <= 0 for slice. var = nanvar(a, axis=axis, dtype=dtype, out=out, ddof=ddof, /Users/nuhahameed/opt/miniconda3/lib/python3.9/site-packages/numpy/lib/nanfunctions.py:1670: RuntimeWarning: Degrees of freedom <= 0 for slice. var = nanvar(a, axis=axis, dtype=dtype, out=out, ddof=ddof, /Users/nuhahameed/opt/miniconda3/lib/python3.9/site-packages/numpy/lib/nanfunctions.py:1670: RuntimeWarning: Degrees of freedom <= 0 for slice. var = nanvar(a, axis=axis, dtype=dtype, out=out, ddof=ddof, /Users/nuhahameed/opt/miniconda3/lib/python3.9/site-packages/numpy/lib/nanfunctions.py:1670: RuntimeWarning: Degrees of freedom <= 0 for slice. var = nanvar(a, axis=axis, dtype=dtype, out=out, ddof=ddof, /Users/nuhahameed/opt/miniconda3/lib/python3.9/site-packages/numpy/lib/nanfunctions.py:1670: RuntimeWarning: Degrees of freedom <= 0 for slice. var = nanvar(a, axis=axis, dtype=dtype, out=out, ddof=ddof, /Users/nuhahameed/opt/miniconda3/lib/python3.9/site-packages/numpy/lib/nanfunctions.py:1670: RuntimeWarning: Degrees of freedom <= 0 for slice. var = nanvar(a, axis=axis, dtype=dtype, out=out, ddof=ddof, /Users/nuhahameed/opt/miniconda3/lib/python3.9/site-packages/numpy/lib/nanfunctions.py:1670: RuntimeWarning: Degrees of freedom <= 0 for slice. var = nanvar(a, axis=axis, dtype=dtype, out=out, ddof=ddof, /Users/nuhahameed/opt/miniconda3/lib/python3.9/site-packages/numpy/lib/nanfunctions.py:1670: RuntimeWarning: Degrees of freedom <= 0 for slice. var = nanvar(a, axis=axis, dtype=dtype, out=out, ddof=ddof, /Users/nuhahameed/opt/miniconda3/lib/python3.9/site-packages/numpy/lib/nanfunctions.py:1670: RuntimeWarning: Degrees of freedom <= 0 for slice. var = nanvar(a, axis=axis, dtype=dtype, out=out, ddof=ddof, /Users/nuhahameed/opt/miniconda3/lib/python3.9/site-packages/numpy/lib/nanfunctions.py:1670: RuntimeWarning: Degrees of freedom <= 0 for slice. var = nanvar(a, axis=axis, dtype=dtype, out=out, ddof=ddof, /Users/nuhahameed/opt/miniconda3/lib/python3.9/site-packages/numpy/lib/nanfunctions.py:1670: RuntimeWarning: Degrees of freedom <= 0 for slice. var = nanvar(a, axis=axis, dtype=dtype, out=out, ddof=ddof, /Users/nuhahameed/opt/miniconda3/lib/python3.9/site-packages/numpy/lib/nanfunctions.py:1670: RuntimeWarning: Degrees of freedom <= 0 for slice. var = nanvar(a, axis=axis, dtype=dtype, out=out, ddof=ddof, /Users/nuhahameed/opt/miniconda3/lib/python3.9/site-packages/numpy/lib/nanfunctions.py:1670: RuntimeWarning: Degrees of freedom <= 0 for slice. var = nanvar(a, axis=axis, dtype=dtype, out=out, ddof=ddof, /Users/nuhahameed/opt/miniconda3/lib/python3.9/site-packages/numpy/lib/nanfunctions.py:1670: RuntimeWarning: Degrees of freedom <= 0 for slice. var = nanvar(a, axis=axis, dtype=dtype, out=out, ddof=ddof, /Users/nuhahameed/opt/miniconda3/lib/python3.9/site-packages/numpy/lib/nanfunctions.py:1670: RuntimeWarning: Degrees of freedom <= 0 for slice. var = nanvar(a, axis=axis, dtype=dtype, out=out, ddof=ddof, /Users/nuhahameed/opt/miniconda3/lib/python3.9/site-packages/numpy/lib/nanfunctions.py:1670: RuntimeWarning: Degrees of freedom <= 0 for slice. var = nanvar(a, axis=axis, dtype=dtype, out=out, ddof=ddof, /Users/nuhahameed/opt/miniconda3/lib/python3.9/site-packages/numpy/lib/nanfunctions.py:1670: RuntimeWarning: Degrees of freedom <= 0 for slice. var = nanvar(a, axis=axis, dtype=dtype, out=out, ddof=ddof, /Users/nuhahameed/opt/miniconda3/lib/python3.9/site-packages/numpy/lib/nanfunctions.py:1670: RuntimeWarning: Degrees of freedom <= 0 for slice. var = nanvar(a, axis=axis, dtype=dtype, out=out, ddof=ddof, /Users/nuhahameed/opt/miniconda3/lib/python3.9/site-packages/numpy/lib/nanfunctions.py:1670: RuntimeWarning: Degrees of freedom <= 0 for slice. var = nanvar(a, axis=axis, dtype=dtype, out=out, ddof=ddof, /Users/nuhahameed/opt/miniconda3/lib/python3.9/site-packages/numpy/lib/nanfunctions.py:1670: RuntimeWarning: Degrees of freedom <= 0 for slice. var = nanvar(a, axis=axis, dtype=dtype, out=out, ddof=ddof, /Users/nuhahameed/opt/miniconda3/lib/python3.9/site-packages/numpy/lib/nanfunctions.py:1670: RuntimeWarning: Degrees of freedom <= 0 for slice. var = nanvar(a, axis=axis, dtype=dtype, out=out, ddof=ddof, /Users/nuhahameed/opt/miniconda3/lib/python3.9/site-packages/numpy/lib/nanfunctions.py:1670: RuntimeWarning: Degrees of freedom <= 0 for slice. var = nanvar(a, axis=axis, dtype=dtype, out=out, ddof=ddof, /Users/nuhahameed/opt/miniconda3/lib/python3.9/site-packages/numpy/lib/nanfunctions.py:1670: RuntimeWarning: Degrees of freedom <= 0 for slice. var = nanvar(a, axis=axis, dtype=dtype, out=out, ddof=ddof, /Users/nuhahameed/opt/miniconda3/lib/python3.9/site-packages/numpy/lib/nanfunctions.py:1670: RuntimeWarning: Degrees of freedom <= 0 for slice. var = nanvar(a, axis=axis, dtype=dtype, out=out, ddof=ddof, /Users/nuhahameed/opt/miniconda3/lib/python3.9/site-packages/numpy/lib/nanfunctions.py:1670: RuntimeWarning: Degrees of freedom <= 0 for slice. var = nanvar(a, axis=axis, dtype=dtype, out=out, ddof=ddof, /Users/nuhahameed/opt/miniconda3/lib/python3.9/site-packages/numpy/lib/nanfunctions.py:1670: RuntimeWarning: Degrees of freedom <= 0 for slice. var = nanvar(a, axis=axis, dtype=dtype, out=out, ddof=ddof, /Users/nuhahameed/opt/miniconda3/lib/python3.9/site-packages/numpy/lib/nanfunctions.py:1670: RuntimeWarning: Degrees of freedom <= 0 for slice. var = nanvar(a, axis=axis, dtype=dtype, out=out, ddof=ddof, /Users/nuhahameed/opt/miniconda3/lib/python3.9/site-packages/numpy/lib/nanfunctions.py:1670: RuntimeWarning: Degrees of freedom <= 0 for slice. var = nanvar(a, axis=axis, dtype=dtype, out=out, ddof=ddof, /Users/nuhahameed/opt/miniconda3/lib/python3.9/site-packages/numpy/lib/nanfunctions.py:1670: RuntimeWarning: Degrees of freedom <= 0 for slice. var = nanvar(a, axis=axis, dtype=dtype, out=out, ddof=ddof, /Users/nuhahameed/opt/miniconda3/lib/python3.9/site-packages/numpy/lib/nanfunctions.py:1670: RuntimeWarning: Degrees of freedom <= 0 for slice. var = nanvar(a, axis=axis, dtype=dtype, out=out, ddof=ddof, /Users/nuhahameed/opt/miniconda3/lib/python3.9/site-packages/numpy/lib/nanfunctions.py:1670: RuntimeWarning: Degrees of freedom <= 0 for slice. var = nanvar(a, axis=axis, dtype=dtype, out=out, ddof=ddof, /Users/nuhahameed/opt/miniconda3/lib/python3.9/site-packages/numpy/lib/nanfunctions.py:1670: RuntimeWarning: Degrees of freedom <= 0 for slice. var = nanvar(a, axis=axis, dtype=dtype, out=out, ddof=ddof, /Users/nuhahameed/opt/miniconda3/lib/python3.9/site-packages/numpy/lib/nanfunctions.py:1670: RuntimeWarning: Degrees of freedom <= 0 for slice. var = nanvar(a, axis=axis, dtype=dtype, out=out, ddof=ddof, /Users/nuhahameed/opt/miniconda3/lib/python3.9/site-packages/numpy/lib/nanfunctions.py:1670: RuntimeWarning: Degrees of freedom <= 0 for slice. var = nanvar(a, axis=axis, dtype=dtype, out=out, ddof=ddof, /Users/nuhahameed/opt/miniconda3/lib/python3.9/site-packages/numpy/lib/nanfunctions.py:1670: RuntimeWarning: Degrees of freedom <= 0 for slice. var = nanvar(a, axis=axis, dtype=dtype, out=out, ddof=ddof, /Users/nuhahameed/opt/miniconda3/lib/python3.9/site-packages/numpy/lib/nanfunctions.py:1670: RuntimeWarning: Degrees of freedom <= 0 for slice. var = nanvar(a, axis=axis, dtype=dtype, out=out, ddof=ddof, /Users/nuhahameed/opt/miniconda3/lib/python3.9/site-packages/numpy/lib/nanfunctions.py:1670: RuntimeWarning: Degrees of freedom <= 0 for slice. var = nanvar(a, axis=axis, dtype=dtype, out=out, ddof=ddof, /Users/nuhahameed/opt/miniconda3/lib/python3.9/site-packages/numpy/lib/nanfunctions.py:1670: RuntimeWarning: Degrees of freedom <= 0 for slice. var = nanvar(a, axis=axis, dtype=dtype, out=out, ddof=ddof, /Users/nuhahameed/opt/miniconda3/lib/python3.9/site-packages/numpy/lib/nanfunctions.py:1670: RuntimeWarning: Degrees of freedom <= 0 for slice. var = nanvar(a, axis=axis, dtype=dtype, out=out, ddof=ddof, /Users/nuhahameed/opt/miniconda3/lib/python3.9/site-packages/numpy/lib/nanfunctions.py:1670: RuntimeWarning: Degrees of freedom <= 0 for slice. var = nanvar(a, axis=axis, dtype=dtype, out=out, ddof=ddof, /Users/nuhahameed/opt/miniconda3/lib/python3.9/site-packages/numpy/lib/nanfunctions.py:1670: RuntimeWarning: Degrees of freedom <= 0 for slice. var = nanvar(a, axis=axis, dtype=dtype, out=out, ddof=ddof, /Users/nuhahameed/opt/miniconda3/lib/python3.9/site-packages/numpy/lib/nanfunctions.py:1670: RuntimeWarning: Degrees of freedom <= 0 for slice. var = nanvar(a, axis=axis, dtype=dtype, out=out, ddof=ddof, /Users/nuhahameed/opt/miniconda3/lib/python3.9/site-packages/numpy/lib/nanfunctions.py:1670: RuntimeWarning: Degrees of freedom <= 0 for slice. var = nanvar(a, axis=axis, dtype=dtype, out=out, ddof=ddof, /Users/nuhahameed/opt/miniconda3/lib/python3.9/site-packages/numpy/lib/nanfunctions.py:1670: RuntimeWarning: Degrees of freedom <= 0 for slice. var = nanvar(a, axis=axis, dtype=dtype, out=out, ddof=ddof, /Users/nuhahameed/opt/miniconda3/lib/python3.9/site-packages/numpy/lib/nanfunctions.py:1670: RuntimeWarning: Degrees of freedom <= 0 for slice. var = nanvar(a, axis=axis, dtype=dtype, out=out, ddof=ddof, /Users/nuhahameed/opt/miniconda3/lib/python3.9/site-packages/numpy/lib/nanfunctions.py:1670: RuntimeWarning: Degrees of freedom <= 0 for slice. var = nanvar(a, axis=axis, dtype=dtype, out=out, ddof=ddof, /Users/nuhahameed/opt/miniconda3/lib/python3.9/site-packages/numpy/lib/nanfunctions.py:1670: RuntimeWarning: Degrees of freedom <= 0 for slice. var = nanvar(a, axis=axis, dtype=dtype, out=out, ddof=ddof, /Users/nuhahameed/opt/miniconda3/lib/python3.9/site-packages/numpy/lib/nanfunctions.py:1670: RuntimeWarning: Degrees of freedom <= 0 for slice. var = nanvar(a, axis=axis, dtype=dtype, out=out, ddof=ddof, /Users/nuhahameed/opt/miniconda3/lib/python3.9/site-packages/numpy/lib/nanfunctions.py:1670: RuntimeWarning: Degrees of freedom <= 0 for slice. var = nanvar(a, axis=axis, dtype=dtype, out=out, ddof=ddof, /Users/nuhahameed/opt/miniconda3/lib/python3.9/site-packages/numpy/lib/nanfunctions.py:1670: RuntimeWarning: Degrees of freedom <= 0 for slice. var = nanvar(a, axis=axis, dtype=dtype, out=out, ddof=ddof, /Users/nuhahameed/opt/miniconda3/lib/python3.9/site-packages/numpy/lib/nanfunctions.py:1670: RuntimeWarning: Degrees of freedom <= 0 for slice. var = nanvar(a, axis=axis, dtype=dtype, out=out, ddof=ddof, /Users/nuhahameed/opt/miniconda3/lib/python3.9/site-packages/numpy/lib/nanfunctions.py:1670: RuntimeWarning: Degrees of freedom <= 0 for slice. var = nanvar(a, axis=axis, dtype=dtype, out=out, ddof=ddof, /Users/nuhahameed/opt/miniconda3/lib/python3.9/site-packages/numpy/lib/nanfunctions.py:1670: RuntimeWarning: Degrees of freedom <= 0 for slice. var = nanvar(a, axis=axis, dtype=dtype, out=out, ddof=ddof, /Users/nuhahameed/opt/miniconda3/lib/python3.9/site-packages/numpy/lib/nanfunctions.py:1670: RuntimeWarning: Degrees of freedom <= 0 for slice. var = nanvar(a, axis=axis, dtype=dtype, out=out, ddof=ddof, /Users/nuhahameed/opt/miniconda3/lib/python3.9/site-packages/numpy/lib/nanfunctions.py:1670: RuntimeWarning: Degrees of freedom <= 0 for slice. var = nanvar(a, axis=axis, dtype=dtype, out=out, ddof=ddof, /Users/nuhahameed/opt/miniconda3/lib/python3.9/site-packages/numpy/lib/nanfunctions.py:1670: RuntimeWarning: Degrees of freedom <= 0 for slice. var = nanvar(a, axis=axis, dtype=dtype, out=out, ddof=ddof, /Users/nuhahameed/opt/miniconda3/lib/python3.9/site-packages/numpy/lib/nanfunctions.py:1670: RuntimeWarning: Degrees of freedom <= 0 for slice. var = nanvar(a, axis=axis, dtype=dtype, out=out, ddof=ddof, /Users/nuhahameed/opt/miniconda3/lib/python3.9/site-packages/numpy/lib/nanfunctions.py:1670: RuntimeWarning: Degrees of freedom <= 0 for slice. var = nanvar(a, axis=axis, dtype=dtype, out=out, ddof=ddof, /Users/nuhahameed/opt/miniconda3/lib/python3.9/site-packages/numpy/lib/nanfunctions.py:1670: RuntimeWarning: Degrees of freedom <= 0 for slice. var = nanvar(a, axis=axis, dtype=dtype, out=out, ddof=ddof, /Users/nuhahameed/opt/miniconda3/lib/python3.9/site-packages/numpy/lib/nanfunctions.py:1670: RuntimeWarning: Degrees of freedom <= 0 for slice. var = nanvar(a, axis=axis, dtype=dtype, out=out, ddof=ddof, /Users/nuhahameed/opt/miniconda3/lib/python3.9/site-packages/numpy/lib/nanfunctions.py:1670: RuntimeWarning: Degrees of freedom <= 0 for slice. var = nanvar(a, axis=axis, dtype=dtype, out=out, ddof=ddof, /Users/nuhahameed/opt/miniconda3/lib/python3.9/site-packages/numpy/lib/nanfunctions.py:1670: RuntimeWarning: Degrees of freedom <= 0 for slice. var = nanvar(a, axis=axis, dtype=dtype, out=out, ddof=ddof, /Users/nuhahameed/opt/miniconda3/lib/python3.9/site-packages/numpy/lib/nanfunctions.py:1670: RuntimeWarning: Degrees of freedom <= 0 for slice. var = nanvar(a, axis=axis, dtype=dtype, out=out, ddof=ddof, /Users/nuhahameed/opt/miniconda3/lib/python3.9/site-packages/numpy/lib/nanfunctions.py:1670: RuntimeWarning: Degrees of freedom <= 0 for slice. var = nanvar(a, axis=axis, dtype=dtype, out=out, ddof=ddof, /Users/nuhahameed/opt/miniconda3/lib/python3.9/site-packages/numpy/lib/nanfunctions.py:1670: RuntimeWarning: Degrees of freedom <= 0 for slice. var = nanvar(a, axis=axis, dtype=dtype, out=out, ddof=ddof, /Users/nuhahameed/opt/miniconda3/lib/python3.9/site-packages/numpy/lib/nanfunctions.py:1670: RuntimeWarning: Degrees of freedom <= 0 for slice. var = nanvar(a, axis=axis, dtype=dtype, out=out, ddof=ddof, /Users/nuhahameed/opt/miniconda3/lib/python3.9/site-packages/numpy/lib/nanfunctions.py:1670: RuntimeWarning: Degrees of freedom <= 0 for slice. var = nanvar(a, axis=axis, dtype=dtype, out=out, ddof=ddof, /Users/nuhahameed/opt/miniconda3/lib/python3.9/site-packages/numpy/lib/nanfunctions.py:1670: RuntimeWarning: Degrees of freedom <= 0 for slice. var = nanvar(a, axis=axis, dtype=dtype, out=out, ddof=ddof, /Users/nuhahameed/opt/miniconda3/lib/python3.9/site-packages/numpy/lib/nanfunctions.py:1670: RuntimeWarning: Degrees of freedom <= 0 for slice. var = nanvar(a, axis=axis, dtype=dtype, out=out, ddof=ddof, /Users/nuhahameed/opt/miniconda3/lib/python3.9/site-packages/numpy/lib/nanfunctions.py:1670: RuntimeWarning: Degrees of freedom <= 0 for slice. var = nanvar(a, axis=axis, dtype=dtype, out=out, ddof=ddof, /Users/nuhahameed/opt/miniconda3/lib/python3.9/site-packages/numpy/lib/nanfunctions.py:1670: RuntimeWarning: Degrees of freedom <= 0 for slice. var = nanvar(a, axis=axis, dtype=dtype, out=out, ddof=ddof, /Users/nuhahameed/opt/miniconda3/lib/python3.9/site-packages/numpy/lib/nanfunctions.py:1670: RuntimeWarning: Degrees of freedom <= 0 for slice. var = nanvar(a, axis=axis, dtype=dtype, out=out, ddof=ddof, /Users/nuhahameed/opt/miniconda3/lib/python3.9/site-packages/numpy/lib/nanfunctions.py:1670: RuntimeWarning: Degrees of freedom <= 0 for slice. var = nanvar(a, axis=axis, dtype=dtype, out=out, ddof=ddof, /Users/nuhahameed/opt/miniconda3/lib/python3.9/site-packages/numpy/lib/nanfunctions.py:1670: RuntimeWarning: Degrees of freedom <= 0 for slice. var = nanvar(a, axis=axis, dtype=dtype, out=out, ddof=ddof,
input_folder = '/Volumes/Seagate Hub/cmip6_downloader/data/Salinity models/'
myfiles = os.listdir(input_folder)
for f in myfiles:
if '._so' in f:
myfiles.remove(f)
def calculate_phi(salinity, depth_range, s0, S=None):
salinity_depth_range = salinity[depth_range]
if S is not None:
salinity_depth_range[0] = S
phi = np.trapz(s0 - salinity_depth_range, depth_range)
return phi
def phi_calculation(file):
ds = xr.open_dataset(input_folder + file)
if 'CESM' in file:
ds['lev'] = ds['lev'] / 100
phi_values = {}
for year in ['1970', '2010']:
ds_mean = ds.so.isel(time=np.arange((int(year) - 1850) * 12, (int(year) + 1 - 1850) * 12, 1)).mean(dim='time')
y = ds_mean.mean(dim=['longitude', 'latitude'])
# Find the nearest depth to 0 and make the max 150
depth_start = np.argmin(np.abs(ds['lev'].values - 0))
depth_end = np.argmin(np.abs(ds['lev'].values - 150))
# Select the depth range
depth_range_values = np.arange(depth_start, depth_end + 1)
# Calculate s0 as the salinity at depth 150
s0 = y.sel(lev=150, method='nearest').values
# Find the nearest depth to 0 and use it to get S
depth_nearest_0 = np.argmin(np.abs(ds['lev'].values - 0))
S = y.sel(lev=ds['lev'].values[depth_nearest_0], method='nearest').values
# Calculate Phi using the updated calculate_phi function
salinity_profile = y.values[depth_range_values]
phi = calculate_phi(salinity_profile, depth_range_values, s0, S)
phi_values[year] = phi
return phi_values
for file in myfiles:
model_name = file.split('so.')[-1].split('.')[0]
phi_values = phi_calculation(file)
results = pd.DataFrame(columns=['Model', 'Phi_1970', 'Phi_2010', 'Phi_ratio'])
for file in myfiles:
model_name = file.split('so.')[-1].split('.')[0]
phi_values = phi_calculation(file)
results = results.append({
'Model': model_name,
'Phi_1970': phi_values['1970'],
'Phi_2010': phi_values['2010'],
'Phi_ratio': phi_values['2010'] / phi_values['1970']
}, ignore_index=True)
results.to_excel('model_phi.xlsx', index=False)
print(results)
/var/folders/ff/60hvmk1s21v54jpfs3v4wp4m0000gn/T/ipykernel_47949/2966426195.py:58: FutureWarning: The frame.append method is deprecated and will be removed from pandas in a future version. Use pandas.concat instead.
results = results.append({
/var/folders/ff/60hvmk1s21v54jpfs3v4wp4m0000gn/T/ipykernel_47949/2966426195.py:58: FutureWarning: The frame.append method is deprecated and will be removed from pandas in a future version. Use pandas.concat instead.
results = results.append({
/var/folders/ff/60hvmk1s21v54jpfs3v4wp4m0000gn/T/ipykernel_47949/2966426195.py:58: FutureWarning: The frame.append method is deprecated and will be removed from pandas in a future version. Use pandas.concat instead.
results = results.append({
/var/folders/ff/60hvmk1s21v54jpfs3v4wp4m0000gn/T/ipykernel_47949/2966426195.py:58: FutureWarning: The frame.append method is deprecated and will be removed from pandas in a future version. Use pandas.concat instead.
results = results.append({
/var/folders/ff/60hvmk1s21v54jpfs3v4wp4m0000gn/T/ipykernel_47949/2966426195.py:58: FutureWarning: The frame.append method is deprecated and will be removed from pandas in a future version. Use pandas.concat instead.
results = results.append({
/var/folders/ff/60hvmk1s21v54jpfs3v4wp4m0000gn/T/ipykernel_47949/2966426195.py:58: FutureWarning: The frame.append method is deprecated and will be removed from pandas in a future version. Use pandas.concat instead.
results = results.append({
/var/folders/ff/60hvmk1s21v54jpfs3v4wp4m0000gn/T/ipykernel_47949/2966426195.py:58: FutureWarning: The frame.append method is deprecated and will be removed from pandas in a future version. Use pandas.concat instead.
results = results.append({
/var/folders/ff/60hvmk1s21v54jpfs3v4wp4m0000gn/T/ipykernel_47949/2966426195.py:58: FutureWarning: The frame.append method is deprecated and will be removed from pandas in a future version. Use pandas.concat instead.
results = results.append({
/var/folders/ff/60hvmk1s21v54jpfs3v4wp4m0000gn/T/ipykernel_47949/2966426195.py:58: FutureWarning: The frame.append method is deprecated and will be removed from pandas in a future version. Use pandas.concat instead.
results = results.append({
/var/folders/ff/60hvmk1s21v54jpfs3v4wp4m0000gn/T/ipykernel_47949/2966426195.py:58: FutureWarning: The frame.append method is deprecated and will be removed from pandas in a future version. Use pandas.concat instead.
results = results.append({
/var/folders/ff/60hvmk1s21v54jpfs3v4wp4m0000gn/T/ipykernel_47949/2966426195.py:58: FutureWarning: The frame.append method is deprecated and will be removed from pandas in a future version. Use pandas.concat instead.
results = results.append({
/var/folders/ff/60hvmk1s21v54jpfs3v4wp4m0000gn/T/ipykernel_47949/2966426195.py:58: FutureWarning: The frame.append method is deprecated and will be removed from pandas in a future version. Use pandas.concat instead.
results = results.append({
/var/folders/ff/60hvmk1s21v54jpfs3v4wp4m0000gn/T/ipykernel_47949/2966426195.py:58: FutureWarning: The frame.append method is deprecated and will be removed from pandas in a future version. Use pandas.concat instead.
results = results.append({
/var/folders/ff/60hvmk1s21v54jpfs3v4wp4m0000gn/T/ipykernel_47949/2966426195.py:58: FutureWarning: The frame.append method is deprecated and will be removed from pandas in a future version. Use pandas.concat instead.
results = results.append({
/var/folders/ff/60hvmk1s21v54jpfs3v4wp4m0000gn/T/ipykernel_47949/2966426195.py:58: FutureWarning: The frame.append method is deprecated and will be removed from pandas in a future version. Use pandas.concat instead.
results = results.append({
/var/folders/ff/60hvmk1s21v54jpfs3v4wp4m0000gn/T/ipykernel_47949/2966426195.py:58: FutureWarning: The frame.append method is deprecated and will be removed from pandas in a future version. Use pandas.concat instead.
results = results.append({
/var/folders/ff/60hvmk1s21v54jpfs3v4wp4m0000gn/T/ipykernel_47949/2966426195.py:58: FutureWarning: The frame.append method is deprecated and will be removed from pandas in a future version. Use pandas.concat instead.
results = results.append({
/var/folders/ff/60hvmk1s21v54jpfs3v4wp4m0000gn/T/ipykernel_47949/2966426195.py:58: FutureWarning: The frame.append method is deprecated and will be removed from pandas in a future version. Use pandas.concat instead.
results = results.append({
/var/folders/ff/60hvmk1s21v54jpfs3v4wp4m0000gn/T/ipykernel_47949/2966426195.py:58: FutureWarning: The frame.append method is deprecated and will be removed from pandas in a future version. Use pandas.concat instead.
results = results.append({
/var/folders/ff/60hvmk1s21v54jpfs3v4wp4m0000gn/T/ipykernel_47949/2966426195.py:58: FutureWarning: The frame.append method is deprecated and will be removed from pandas in a future version. Use pandas.concat instead.
results = results.append({
/var/folders/ff/60hvmk1s21v54jpfs3v4wp4m0000gn/T/ipykernel_47949/2966426195.py:58: FutureWarning: The frame.append method is deprecated and will be removed from pandas in a future version. Use pandas.concat instead.
results = results.append({
/var/folders/ff/60hvmk1s21v54jpfs3v4wp4m0000gn/T/ipykernel_47949/2966426195.py:58: FutureWarning: The frame.append method is deprecated and will be removed from pandas in a future version. Use pandas.concat instead.
results = results.append({
/var/folders/ff/60hvmk1s21v54jpfs3v4wp4m0000gn/T/ipykernel_47949/2966426195.py:58: FutureWarning: The frame.append method is deprecated and will be removed from pandas in a future version. Use pandas.concat instead.
results = results.append({
/var/folders/ff/60hvmk1s21v54jpfs3v4wp4m0000gn/T/ipykernel_47949/2966426195.py:58: FutureWarning: The frame.append method is deprecated and will be removed from pandas in a future version. Use pandas.concat instead.
results = results.append({
/var/folders/ff/60hvmk1s21v54jpfs3v4wp4m0000gn/T/ipykernel_47949/2966426195.py:58: FutureWarning: The frame.append method is deprecated and will be removed from pandas in a future version. Use pandas.concat instead.
results = results.append({
/var/folders/ff/60hvmk1s21v54jpfs3v4wp4m0000gn/T/ipykernel_47949/2966426195.py:58: FutureWarning: The frame.append method is deprecated and will be removed from pandas in a future version. Use pandas.concat instead.
results = results.append({
/var/folders/ff/60hvmk1s21v54jpfs3v4wp4m0000gn/T/ipykernel_47949/2966426195.py:58: FutureWarning: The frame.append method is deprecated and will be removed from pandas in a future version. Use pandas.concat instead.
results = results.append({
/var/folders/ff/60hvmk1s21v54jpfs3v4wp4m0000gn/T/ipykernel_47949/2966426195.py:58: FutureWarning: The frame.append method is deprecated and will be removed from pandas in a future version. Use pandas.concat instead.
results = results.append({
/var/folders/ff/60hvmk1s21v54jpfs3v4wp4m0000gn/T/ipykernel_47949/2966426195.py:58: FutureWarning: The frame.append method is deprecated and will be removed from pandas in a future version. Use pandas.concat instead.
results = results.append({
Model Phi_1970 Phi_2010 Phi_ratio 0 BCC-ESM1 20.710111 20.365960 0.983382 1 CNRM-CM6-1-HR 22.106651 23.435199 1.060097 2 CAMS-CSM1-0 19.176156 19.075460 0.994749 3 GISS-E2-1-G 9.707650 11.073850 1.140734 4 CanESM5 7.367585 15.315710 2.078796 5 CNRM-CM6-1 38.011110 43.646579 1.148258 6 CESM2-WACCM-FV2 12.967559 15.750220 1.214586 7 TaiESM1 10.439663 16.532539 1.583628 8 CESM2 10.933884 15.796343 1.444715 9 CNRM-ESM2-1 39.249865 41.851585 1.066286 10 CESM2-WACCM 11.830647 13.976324 1.181366 11 MPI-ESM1-2-HR 5.881024 9.908315 1.684794 12 BCC-CSM2-MR 19.165132 20.522662 1.070833 13 CESM2-FV2 11.645007 16.831038 1.445344 14 EC-Earth3-CC 16.198829 47.973336 2.961531 15 EC-Earth3-Veg-LR 14.044454 29.706441 2.115172 16 EC-Earth3 40.677680 30.289046 0.744611 17 E3SM-1-0 38.377339 37.991172 0.989938 18 HadGEM3-GC31-LL 21.725981 36.421713 1.676413 19 CMCC-CM2-SR5 17.999363 16.291018 0.905089 20 CMCC-ESM2 15.582287 16.279681 1.044756 21 EC-Earth3-AerChem 11.228786 25.333605 2.256130 22 EC-Earth3-Veg 14.320418 33.341941 2.328280 23 FIO-ESM-2-0 8.481615 14.404881 1.698365 24 GFDL-CM4 3.128994 5.291096 1.690989 25 GFDL-ESM4 3.410776 4.295567 1.259410 26 GISS-E2-2-G 8.499157 7.432849 0.874540 27 MPI-ESM-1-2-HAM 9.447742 9.488862 1.004352 28 MPI-ESM1-2-LR 9.158100 11.342212 1.238490 29 MRI-ESM2-0 9.053513 11.182763 1.235185
/var/folders/ff/60hvmk1s21v54jpfs3v4wp4m0000gn/T/ipykernel_47949/2966426195.py:58: FutureWarning: The frame.append method is deprecated and will be removed from pandas in a future version. Use pandas.concat instead.
results = results.append({
input_folder = '/Volumes/Seagate Hub/cmip6_downloader/data/Salinity models/'
myfiles = os.listdir(input_folder)
for f in myfiles:
if '._so' in f:
myfiles.remove(f)
def calculate_phi(salinity, depth_range):
salinity_depth_range = salinity[depth_range]
Phi = np.trapz(salinity_depth_range - 33)
return Phi
def phi_calculation(file):
ds = xr.open_dataset(input_folder + file)
if 'CESM' in file:
ds['lev'] = ds['lev'] / 100
phi_values = {}
for year in ['1970', '2010']:
ds_mean = ds.so.isel(time=np.arange((int(year) - 1850) * 12, (int(year) + 1 - 1850) * 12, 1)).mean(dim='time')
y = ds_mean.mean(dim=['longitude', 'latitude'])
# Find the nearest depth to 0 and make the max 150
depth_start = np.argmin(np.abs(ds['lev'].values - 0))
depth_end = np.argmin(np.abs(ds['lev'].values - 150))
# Select the depth range
depth_range_values = np.arange(depth_start, depth_end + 1)
salinity_profile = y.values[depth_range_values]
phi = calculate_phi(salinity_profile, depth_range_values)
phi_values[year] = phi
return phi_values
depth_range = np.arange(0, 151)
for file in myfiles:
model_name = file.split('so.')[-1].split('.')[0]
phi_values = phi_calculation(file)
results = pd.DataFrame(columns=['Model', 'Phi_1970', 'Phi_2010', 'Phi_ratio'])
for file in myfiles:
model_name = file.split('so.')[-1].split('.')[0]
phi_values = phi_calculation(file)
results = results.append({
'Model': model_name,
'Phi_1970': phi_values['1970'],
'Phi_2010': phi_values['2010'],
'Phi_ratio': phi_values['2010'] / phi_values['1970']
}, ignore_index=True)
results.to_excel('model_phi.xlsx', index=False)
print(results)
/var/folders/ff/60hvmk1s21v54jpfs3v4wp4m0000gn/T/ipykernel_47949/113800120.py:48: FutureWarning: The frame.append method is deprecated and will be removed from pandas in a future version. Use pandas.concat instead.
results = results.append({
/var/folders/ff/60hvmk1s21v54jpfs3v4wp4m0000gn/T/ipykernel_47949/113800120.py:48: FutureWarning: The frame.append method is deprecated and will be removed from pandas in a future version. Use pandas.concat instead.
results = results.append({
/var/folders/ff/60hvmk1s21v54jpfs3v4wp4m0000gn/T/ipykernel_47949/113800120.py:48: FutureWarning: The frame.append method is deprecated and will be removed from pandas in a future version. Use pandas.concat instead.
results = results.append({
/var/folders/ff/60hvmk1s21v54jpfs3v4wp4m0000gn/T/ipykernel_47949/113800120.py:48: FutureWarning: The frame.append method is deprecated and will be removed from pandas in a future version. Use pandas.concat instead.
results = results.append({
/var/folders/ff/60hvmk1s21v54jpfs3v4wp4m0000gn/T/ipykernel_47949/113800120.py:48: FutureWarning: The frame.append method is deprecated and will be removed from pandas in a future version. Use pandas.concat instead.
results = results.append({
/var/folders/ff/60hvmk1s21v54jpfs3v4wp4m0000gn/T/ipykernel_47949/113800120.py:48: FutureWarning: The frame.append method is deprecated and will be removed from pandas in a future version. Use pandas.concat instead.
results = results.append({
/var/folders/ff/60hvmk1s21v54jpfs3v4wp4m0000gn/T/ipykernel_47949/113800120.py:48: FutureWarning: The frame.append method is deprecated and will be removed from pandas in a future version. Use pandas.concat instead.
results = results.append({
/var/folders/ff/60hvmk1s21v54jpfs3v4wp4m0000gn/T/ipykernel_47949/113800120.py:48: FutureWarning: The frame.append method is deprecated and will be removed from pandas in a future version. Use pandas.concat instead.
results = results.append({
/var/folders/ff/60hvmk1s21v54jpfs3v4wp4m0000gn/T/ipykernel_47949/113800120.py:48: FutureWarning: The frame.append method is deprecated and will be removed from pandas in a future version. Use pandas.concat instead.
results = results.append({
/var/folders/ff/60hvmk1s21v54jpfs3v4wp4m0000gn/T/ipykernel_47949/113800120.py:48: FutureWarning: The frame.append method is deprecated and will be removed from pandas in a future version. Use pandas.concat instead.
results = results.append({
/var/folders/ff/60hvmk1s21v54jpfs3v4wp4m0000gn/T/ipykernel_47949/113800120.py:48: FutureWarning: The frame.append method is deprecated and will be removed from pandas in a future version. Use pandas.concat instead.
results = results.append({
/var/folders/ff/60hvmk1s21v54jpfs3v4wp4m0000gn/T/ipykernel_47949/113800120.py:48: FutureWarning: The frame.append method is deprecated and will be removed from pandas in a future version. Use pandas.concat instead.
results = results.append({
/var/folders/ff/60hvmk1s21v54jpfs3v4wp4m0000gn/T/ipykernel_47949/113800120.py:48: FutureWarning: The frame.append method is deprecated and will be removed from pandas in a future version. Use pandas.concat instead.
results = results.append({
/var/folders/ff/60hvmk1s21v54jpfs3v4wp4m0000gn/T/ipykernel_47949/113800120.py:48: FutureWarning: The frame.append method is deprecated and will be removed from pandas in a future version. Use pandas.concat instead.
results = results.append({
/var/folders/ff/60hvmk1s21v54jpfs3v4wp4m0000gn/T/ipykernel_47949/113800120.py:48: FutureWarning: The frame.append method is deprecated and will be removed from pandas in a future version. Use pandas.concat instead.
results = results.append({
/var/folders/ff/60hvmk1s21v54jpfs3v4wp4m0000gn/T/ipykernel_47949/113800120.py:48: FutureWarning: The frame.append method is deprecated and will be removed from pandas in a future version. Use pandas.concat instead.
results = results.append({
/var/folders/ff/60hvmk1s21v54jpfs3v4wp4m0000gn/T/ipykernel_47949/113800120.py:48: FutureWarning: The frame.append method is deprecated and will be removed from pandas in a future version. Use pandas.concat instead.
results = results.append({
/var/folders/ff/60hvmk1s21v54jpfs3v4wp4m0000gn/T/ipykernel_47949/113800120.py:48: FutureWarning: The frame.append method is deprecated and will be removed from pandas in a future version. Use pandas.concat instead.
results = results.append({
/var/folders/ff/60hvmk1s21v54jpfs3v4wp4m0000gn/T/ipykernel_47949/113800120.py:48: FutureWarning: The frame.append method is deprecated and will be removed from pandas in a future version. Use pandas.concat instead.
results = results.append({
/var/folders/ff/60hvmk1s21v54jpfs3v4wp4m0000gn/T/ipykernel_47949/113800120.py:48: FutureWarning: The frame.append method is deprecated and will be removed from pandas in a future version. Use pandas.concat instead.
results = results.append({
/var/folders/ff/60hvmk1s21v54jpfs3v4wp4m0000gn/T/ipykernel_47949/113800120.py:48: FutureWarning: The frame.append method is deprecated and will be removed from pandas in a future version. Use pandas.concat instead.
results = results.append({
/var/folders/ff/60hvmk1s21v54jpfs3v4wp4m0000gn/T/ipykernel_47949/113800120.py:48: FutureWarning: The frame.append method is deprecated and will be removed from pandas in a future version. Use pandas.concat instead.
results = results.append({
/var/folders/ff/60hvmk1s21v54jpfs3v4wp4m0000gn/T/ipykernel_47949/113800120.py:48: FutureWarning: The frame.append method is deprecated and will be removed from pandas in a future version. Use pandas.concat instead.
results = results.append({
/var/folders/ff/60hvmk1s21v54jpfs3v4wp4m0000gn/T/ipykernel_47949/113800120.py:48: FutureWarning: The frame.append method is deprecated and will be removed from pandas in a future version. Use pandas.concat instead.
results = results.append({
/var/folders/ff/60hvmk1s21v54jpfs3v4wp4m0000gn/T/ipykernel_47949/113800120.py:48: FutureWarning: The frame.append method is deprecated and will be removed from pandas in a future version. Use pandas.concat instead.
results = results.append({
/var/folders/ff/60hvmk1s21v54jpfs3v4wp4m0000gn/T/ipykernel_47949/113800120.py:48: FutureWarning: The frame.append method is deprecated and will be removed from pandas in a future version. Use pandas.concat instead.
results = results.append({
Model Phi_1970 Phi_2010 Phi_ratio 0 BCC-ESM1 -25.370138 -25.460762 1.003572 1 CNRM-CM6-1-HR -22.007671 -23.671782 1.075615 2 CAMS-CSM1-0 -34.605343 -34.153458 0.986942 3 GISS-E2-1-G -7.055215 -8.449552 1.197632 4 CanESM5 -6.361787 -16.368155 2.572886 5 CNRM-CM6-1 -36.197334 -51.601128 1.425550 6 CESM2-WACCM-FV2 -14.957618 -20.495430 1.370234 7 TaiESM1 -6.863827 -16.935165 2.467307 8 CESM2 -11.750298 -19.035990 1.620043 9 CNRM-ESM2-1 -45.352631 -54.328934 1.197922 10 CESM2-WACCM -9.820188 -15.971029 1.626347 11 MPI-ESM1-2-HR -3.259603 -11.067601 3.395384 12 BCC-CSM2-MR -27.402544 -28.595207 1.043524 13 CESM2-FV2 -9.887156 -19.254217 1.947397 14 EC-Earth3-CC -21.264727 -59.704468 2.807676 15 EC-Earth3-Veg-LR -5.669782 -25.642700 4.522696 16 EC-Earth3 -36.740936 -42.166641 1.147675 17 E3SM-1-0 -32.565140 -39.827255 1.223003 18 HadGEM3-GC31-LL -21.562834 -35.387524 1.641135 19 CMCC-CM2-SR5 -30.648678 -36.354385 1.186165 20 CMCC-ESM2 -28.809265 -33.947563 1.178356 21 EC-Earth3-AerChem -8.978821 -27.347092 3.045733 22 EC-Earth3-Veg -12.333816 -38.470356 3.119096 23 FIO-ESM-2-0 -6.334970 -15.749584 2.486134 24 GFDL-CM4 2.249943 -1.767567 -0.785605 25 GFDL-ESM4 -0.138193 -0.557011 4.030668 26 GISS-E2-2-G -2.429974 -2.009784 0.827080 27 MPI-ESM-1-2-HAM -3.890644 -7.132752 1.833309 28 MPI-ESM1-2-LR -7.579884 -12.229815 1.613457 29 MRI-ESM2-0 -5.015831 -11.696031 2.331823
/var/folders/ff/60hvmk1s21v54jpfs3v4wp4m0000gn/T/ipykernel_47949/113800120.py:48: FutureWarning: The frame.append method is deprecated and will be removed from pandas in a future version. Use pandas.concat instead.
results = results.append({
/var/folders/ff/60hvmk1s21v54jpfs3v4wp4m0000gn/T/ipykernel_47949/113800120.py:48: FutureWarning: The frame.append method is deprecated and will be removed from pandas in a future version. Use pandas.concat instead.
results = results.append({
/var/folders/ff/60hvmk1s21v54jpfs3v4wp4m0000gn/T/ipykernel_47949/113800120.py:48: FutureWarning: The frame.append method is deprecated and will be removed from pandas in a future version. Use pandas.concat instead.
results = results.append({
/var/folders/ff/60hvmk1s21v54jpfs3v4wp4m0000gn/T/ipykernel_47949/113800120.py:48: FutureWarning: The frame.append method is deprecated and will be removed from pandas in a future version. Use pandas.concat instead.
results = results.append({
input_folder = '/Volumes/Seagate Hub/cmip6_downloader/data/Salinity models/'
myfiles = os.listdir(input_folder)
for f in myfiles:
if '._so' in f:
myfiles.remove(f)
colors = sns.color_palette("husl", len(myfiles))
linestyles = ['-', '--', '-.', ':'] * (len(myfiles) // 4 + 1)
fig, axs = plt.subplots(nrows=1, ncols=2, figsize=(29, 12), sharex=True)
plt.subplots_adjust(right=0.7, hspace=0.3)
def plot_salinity_profile(ax, file, color, linestyle, start_year, end_year):
ds = xr.open_dataset(input_folder + file)
if 'CESM' in file:
ds['lev'] = ds['lev'] / 100
model_name = file.split('so.')[-1].split('.')[0]
#start_year = 2000
#end_year = 2010
ds_mean = ds.so.isel(time=np.arange((start_year-1850)*12, (end_year+1-1850)*12, 1)).mean(dim='time')
y = ds_mean.mean(dim=['longitude', 'latitude'])
ax.plot(y, ds['lev'], label=model_name, color=color, linewidth=2, linestyle=linestyle)
for file, color, linestyle in zip(myfiles, colors, linestyles):
plot_salinity_profile(axs[0], file, color, linestyle, 1970, 1980)
plot_salinity_profile(axs[1], file, color, linestyle, 2000, 2010)
axs[0].plot(AJX_avg_salinity, depth_values, label='1975 AIDJEX', color='red', linewidth = 4)
axs[1].plot(ITP_avg_salinity, depth_values, label='2006-2012 ITP', color='orange', linewidth=4)
for i, ax in enumerate(axs):
ax.set_ylabel('Depth (m)', fontsize=14)
ax.set_xlabel('Salinity (g/kg)', fontsize=14)
ax.set_ylim(0, 150)
ax.set_xlim(27, 35)
ax.tick_params(axis='y', labelsize=12)
ax.tick_params(axis='x', labelsize=12)
ax.invert_yaxis()
# Add subplot labels
label = chr(97 + i) # Convert index to ASCII character ('a', 'b', etc.)
ax.text(-0.1, 1, label, fontsize=19, fontweight='bold', transform=ax.transAxes)
axs[0].set_title('1970-1980 Salinity Profiles', fontsize=25)
axs[1].set_title('2000-2010 Salinity Profiles', fontsize=25)
# Get the handles and labels from the first subplot
handles, labels = axs[0].get_legend_handles_labels()
# Add the handle and label for the ITP line from the second subplot
handles_itp, labels_itp = axs[1].get_legend_handles_labels()
handles += handles_itp[-1:]
labels += labels_itp[-1:]
# Create a shared legend for all subplots
fig.legend(handles, labels, fontsize=15, bbox_to_anchor=(0.68, -0.1, 0.15, 1), ncol=1, borderaxespad=0)
plt.savefig('multi_model_Salinity_Profiles_13_04_.png', bbox_inches='tight')
# Define input folder path
input_folder = '/Volumes/Seagate Hub/cmip6_downloader/data/Salinity models/'
# Create list of all files in folder
all_files = os.listdir(input_folder)
# Remove any hidden files from the list
all_files = [f for f in all_files if not f.startswith('.')]
def calc_mean_so(ds, year):
#Year = ['1970', '1980', '1990', '2000', '2010']
salinity = []
#for j, year in enumerate(Year):
ds_mean = ds.so.isel(time=np.arange((int(year)-1850)*12, (int(year)+1-1850)*12,1)).mean(dim='time')
y = ds_mean.mean(dim=['longitude', 'latitude']) #works out the average within the longitudanal range
salinity.append(y)
#print(len(y))
#else:
# print(f"No data found for {year}")
return salinity
# Loop over all files and calculate decadal average salinity data
all_salinity_depth_data = []
#all_salinity_lat_lon_data = []
all_y = []
for file in all_files:
salinity_depth_data = calc_mean_so(xr.open_dataset(input_folder+file), '1980')
all_salinity_depth_data.append(salinity_depth_data)
# all_salinity_lat_lon_data.append(salinity_lat_lon_data)
#all_salinity_depth_data
new_list = []
def interp(array):
x_out = np.arange(0, 150, 1)
y_in = array.values
x_in = array.lev.values
f = interpolate.interp1d(x_in, y_in, fill_value='extrapolate') #option: change to fill_value='extrapolate'
y_out = f(x_out)
return y_out
for i in range(len(all_salinity_depth_data)):
array = all_salinity_depth_data[i]
new_array = interp(array[0])
new_list.append(new_array)
mean_out_1980 = np.mean(new_list, axis=0)
# Loop over all files and calculate decadal average salinity data
all_salinity_depth_data = []
#all_salinity_lat_lon_data = []
all_y = []
for file in all_files:
salinity_depth_data = calc_mean_so(xr.open_dataset(input_folder+file), '1990')
all_salinity_depth_data.append(salinity_depth_data)
# all_salinity_lat_lon_data.append(salinity_lat_lon_data)
#all_salinity_depth_data
new_list = []
def interp(array):
x_out = np.arange(0, 150, 1)
y_in = array.values
x_in = array.lev.values
f = interpolate.interp1d(x_in, y_in, fill_value='extrapolate') #option: change to fill_value='extrapolate'
y_out = f(x_out)
return y_out
for i in range(len(all_salinity_depth_data)):
array = all_salinity_depth_data[i]
new_array = interp(array[0])
new_list.append(new_array)
mean_out_1990 = np.mean(new_list, axis=0)
# Loop over all files and calculate decadal average salinity data
all_salinity_depth_data = []
#all_salinity_lat_lon_data = []
all_y = []
for file in all_files:
salinity_depth_data = calc_mean_so(xr.open_dataset(input_folder+file), '2000')
all_salinity_depth_data.append(salinity_depth_data)
# all_salinity_lat_lon_data.append(salinity_lat_lon_data)
#all_salinity_depth_data
new_list = []
def interp(array):
x_out = np.arange(0, 150, 1)
y_in = array.values
x_in = array.lev.values
f = interpolate.interp1d(x_in, y_in, fill_value='extrapolate') #option: change to fill_value='extrapolate'
y_out = f(x_out)
return y_out
for i in range(len(all_salinity_depth_data)):
array = all_salinity_depth_data[i]
new_array = interp(array[0])
new_list.append(new_array)
mean_out_2000 = np.mean(new_list, axis=0)
# Loop over all files and calculate decadal average salinity data
all_salinity_depth_data = []
#all_salinity_lat_lon_data = []
all_y = []
for file in all_files:
salinity_depth_data = calc_mean_so(xr.open_dataset(input_folder+file), '2010')
all_salinity_depth_data.append(salinity_depth_data)
# all_salinity_lat_lon_data.append(salinity_lat_lon_data)
#all_salinity_depth_data
new_list = []
def interp(array):
x_out = np.arange(0, 150, 1)
y_in = array.values
x_in = array.lev.values
f = interpolate.interp1d(x_in, y_in, fill_value='extrapolate') #option: change to fill_value='extrapolate'
y_out = f(x_out)
return y_out
for i in range(len(all_salinity_depth_data)):
array = all_salinity_depth_data[i]
new_array = interp(array[0])
new_list.append(new_array)
mean_out_2010 = np.mean(new_list, axis=0)
def plot_salinity_depth(years, salinity_depth_data):
for i, data in enumerate(salinity_depth_data):
Year = ['1970', '1980', '1990', '2000', '2010']
plt.plot(data, np.arange(0, 150, 1), label=years[i], color=plt.cm.viridis(i/len(Year)), linewidth = 4)
plt.ylim(150, 0)
plt.xlim(27,35)
plt.xlabel('Salinity (g/kg)', fontsize = 15)
plt.ylabel('Depth (m)', fontsize = 15)
plt.title(f"Multi Model Mean", fontsize = 20)
plt.plot(AJX_avg_salinity, depth_values, label='1975 AIDJEX', color='red', linewidth = 4)
plt.plot(ITP_avg_salinity, depth_values, label='2006-2012 ITP', color='orange', linewidth = 4)
plt.legend()
plt.savefig('Multi Model Mean.png')
plt.show()
# Loop over all files and calculate decadal average salinity data
all_salinity_depth_data = {}
for year in ['1970', '1980', '1990', '2000', '2010']:
salinity_depth_data = []
for file in all_files:
salinity_depth_data.append(calc_mean_so(xr.open_dataset(input_folder+file), year))
all_salinity_depth_data[year] = salinity_depth_data
# Interpolate the data
interpolated_salinity_data = {}
for year, salinity_depth_data in all_salinity_depth_data.items():
interpolated_salinity_data[year] = [interp(array[0]) for array in salinity_depth_data]
# Calculate the mean salinity data for each decade
mean_salinity_data = {}
for year, interpolated_salinity in interpolated_salinity_data.items():
mean_salinity_data[year] = np.mean(interpolated_salinity, axis=0)
# Plot the mean salinity data for all decades
plot_salinity_depth(['1970', '1980', '1990', '2000', '2010'], list(mean_salinity_data.values()))
# Define a function to calculate maximum salinity gradient and depth for a given model and time period
def max_gradient_depth(file, year='2010', max_depth=150):
ds = xr.open_dataset(os.path.join(input_folder, file))
if 'CESM' in file:
ds['lev'] = ds['lev']/100
model_name = file.split('so.')[-1].split('.')[0]
# Extract salinity data for the specified year
ds_mean = ds.so.sel(lev=slice(0, max_depth)).isel(time=np.arange((int(year)-1850)*12, (int(year)+1-1850)*12,1)).mean(dim='time')
y = ds_mean.mean(dim=['longitude', 'latitude'])
# Calculate salinity gradient (derivative)
d_values = ds.lev.sel(lev=slice(0, max_depth)).values
ds_derivative = np.gradient(y, d_values, axis=0)
# Find maximum salinity gradient and depth
max_gradient = ds_derivative[:-2].max()
depth_max_gradient = d_values[ds_derivative[:-2].argmax()]
# Extract salinity value at depth 0 for the specified year
salinity_at_depth_0 = ds['so'].sel(time='2010', lev=ds['lev'].sel(method='nearest', lev=0)).mean().values
return {'Model': model_name, 'Max Gradient': max_gradient, 'Mixed layer depth': depth_max_gradient, 'Surface Salinity': salinity_at_depth_0}
# Create an empty list to store results
results = []
# Loop through each file and calculate maximum salinity gradient and depth for 2010
for file in myfiles:
result = max_gradient_depth(file)
results.append(result)
# Convert results to a pandas DataFrame
df = pd.DataFrame(results)
# Sort results by maximum gradient in descending order
#df = df.sort_values('Max Gradient', ascending=False)
# Print table
print(df[['Model', 'Surface Salinity', 'Max Gradient', 'Mixed layer depth']])
surface_salinity_1970 = [data[0] for data in interpolated_salinity_data['1970']] # Salinity data at depth=0m for 1975
surface_salinity_2010 = [data[0] for data in interpolated_salinity_data['2010']] # Salinity data at depth=0m for 2010
surface_salinity_diff = np.array(surface_salinity_1970) - np.array(surface_salinity_2010)
surface_salinity_diff[0]
input_folder = '/Volumes/Seagate Hub/cmip6_downloader/data/Salinity models/'
myfiles = os.listdir(input_folder)
# Remove any filenames that contain '._so'
myfiles = [f for f in myfiles if '._so' not in f]
# Create a list to store the differences in average salinity at depth 0
salinity_diffs = []
# Iterate over each file
for file in myfiles:
ds = xr.open_dataset(input_folder+file)
# Extract the model name from the filename
model_name = file.split('so.')[-1].split('.')[0]
# Calculate the mean salinity at depth 0 for each year
salinity_1970 = ds['so'].sel(time='1970', lev=ds['lev'].sel(method='nearest', lev=0)).mean().values
salinity_2010 = ds['so'].sel(time='2010', lev=ds['lev'].sel(method='nearest', lev=0)).mean().values
# Calculate the difference in average salinity at depth 0 between 1970 and 2010
salinity_diff = salinity_2010 - salinity_1970
# Add the model name and salinity difference to the list
salinity_diffs.append([model_name, salinity_diff])
# Create a pandas DataFrame from the list
df = pd.DataFrame(salinity_diffs, columns=['Model name', 'Salinity difference'])
df = df.sort_values('Salinity difference', ascending=False)
df.to_excel('So_differences.xlsx', index=False)
print(df)
Model name Salinity difference 16 EC-Earth3 0.319889 0 BCC-ESM1 0.068560 2 CAMS-CSM1-0 0.038019 3 GISS-E2-1-G -0.058065 26 GISS-E2-2-G -0.059908 25 GFDL-ESM4 -0.107124 1 CNRM-CM6-1-HR -0.109035 19 CMCC-CM2-SR5 -0.175518 12 BCC-CSM2-MR -0.190695 27 MPI-ESM-1-2-HAM -0.193058 20 CMCC-ESM2 -0.204292 9 CNRM-ESM2-1 -0.367008 28 MPI-ESM1-2-LR -0.397720 29 MRI-ESM2-0 -0.400244 6 CESM2-WACCM-FV2 -0.457796 5 CNRM-CM6-1 -0.463995 10 CESM2-WACCM -0.520042 11 MPI-ESM1-2-HR -0.575769 8 CESM2 -0.642033 24 GFDL-CM4 -0.655903 17 E3SM-1-0 -0.715244 13 CESM2-FV2 -0.731340 18 HadGEM3-GC31-LL -0.733768 15 EC-Earth3-Veg-LR -0.777430 23 FIO-ESM-2-0 -0.934566 7 TaiESM1 -1.015944 4 CanESM5 -1.130690 22 EC-Earth3-Veg -1.234190 21 EC-Earth3-AerChem -1.273705 14 EC-Earth3-CC -2.438274
--------------------------------------------------------------------------- KeyError Traceback (most recent call last) File ~/opt/miniconda3/lib/python3.9/site-packages/pandas/core/indexes/base.py:3629, in Index.get_loc(self, key, method, tolerance) 3628 try: -> 3629 return self._engine.get_loc(casted_key) 3630 except KeyError as err: File ~/opt/miniconda3/lib/python3.9/site-packages/pandas/_libs/index.pyx:136, in pandas._libs.index.IndexEngine.get_loc() File ~/opt/miniconda3/lib/python3.9/site-packages/pandas/_libs/index.pyx:163, in pandas._libs.index.IndexEngine.get_loc() File pandas/_libs/hashtable_class_helper.pxi:5198, in pandas._libs.hashtable.PyObjectHashTable.get_item() File pandas/_libs/hashtable_class_helper.pxi:5206, in pandas._libs.hashtable.PyObjectHashTable.get_item() KeyError: 'Difference in salinity between 2010 and 1970 (g/kg)' The above exception was the direct cause of the following exception: KeyError Traceback (most recent call last) Cell In [4], line 37 34 print(df) 36 # Plot a histogram of the salinity difference values ---> 37 plt.hist(df['Difference in salinity between 2010 and 1970 (g/kg)'], bins=10) 38 plt.title('Salinity Difference Histogram') 39 plt.xlabel('Salinity Difference') File ~/opt/miniconda3/lib/python3.9/site-packages/pandas/core/frame.py:3505, in DataFrame.__getitem__(self, key) 3503 if self.columns.nlevels > 1: 3504 return self._getitem_multilevel(key) -> 3505 indexer = self.columns.get_loc(key) 3506 if is_integer(indexer): 3507 indexer = [indexer] File ~/opt/miniconda3/lib/python3.9/site-packages/pandas/core/indexes/base.py:3631, in Index.get_loc(self, key, method, tolerance) 3629 return self._engine.get_loc(casted_key) 3630 except KeyError as err: -> 3631 raise KeyError(key) from err 3632 except TypeError: 3633 # If we have a listlike key, _check_indexing_error will raise 3634 # InvalidIndexError. Otherwise we fall through and re-raise 3635 # the TypeError. 3636 self._check_indexing_error(key) KeyError: 'Difference in salinity between 2010 and 1970 (g/kg)'
# Plot a histogram of the salinity difference values
plt.hist(df['Salinity difference'], bins=10)
plt.title('Surface Salinity Difference', fontsize=15)
plt.xlabel('Difference in surface salinity between 2010 and 1970 (g/kg)', fontsize=12)
plt.ylabel('Frequency', fontsize=12)
plt.savefig('hist_sal_diff')
obs_so_diff = ITP_avg_salinity[0:1] - AJX_avg_salinity[0:1]
#obs_so_diff[0:1]
obs_so_diff
array([-3.42210326])
excel_file_path = '/Volumes/Seagate Hub/ocean_components.xlsx'
model_info_df = pd.read_excel(excel_file_path, engine='openpyxl', sheet_name='Sheet1', header=None, names=['Model Name', 'Ocean Component'])
# Display the DataFrame
print(model_info_df)
Model Name Ocean Component 0 BCC-ESM1 MOM4-L40 1 CNRM-CM6-1-HR NEMO3.6 2 CAMS-CSM1-0 MOM4 3 GISS-E2-1-G GISS Ocean 4 CanESM5 Nemo3.4.1 5 CNRM-CM6-1 NEMO3.6 6 CESM2-WACCM-FV2 POP2 7 MCM-UA-1-0 NaN 8 CESM2 POP2 9 CNRM-ESM2-1 NEMO3.6 10 CESM2-WACCM POP2 11 MPI-ESM1-2-HR MPIOM1.6.3 12 BCC-CSM2-MR MOM4-L40 13 CESM2-FV2 POP2 14 EC-Earth3-CC NEMO3.6 15 EC-Earth3-Veg-LR NEMO3.6 16 EC-Earth3 NEMO3.6 17 E3SM-1-0 MPAS-Ocean 18 HadGEM3-GC31-LL NEMO-HadGEM3-GO6.0 19 CMCC-CM2-SR5 NEMO3.6 20 CMCC-ESM2 NEMO3.6 21 EC-Earth3-AerChem NEMO-LIM3.3.6 22 EC-Earth3-Veg NEMO3.6 23 FIO-ESM-2-0 POP2 24 GFDL-CM4 MOM6 25 GFDL-ESM4 MOM6 26 GISS-E2-2-G GISS Ocean 27 MPI-ESM-1-2-HAM MPIOM1.6.3 28 MPI-ESM1-2-LR MPIOM1.6.3 29 MRI-ESM2-0 MRI.COM4.4 30 TaiESM1 POP2
# Read the Excel sheet into a pandas DataFrame
excel_file_path = '/Volumes/Seagate Hub/ocean_components.xlsx'
model_info_df = pd.read_excel(excel_file_path, engine='openpyxl', sheet_name='Sheet1', header=None, names=['Model Name', 'Ocean Component'])
ocean_component_colors = {
'MOM4-L40':'magenta',
'MOM4':'navy',
'MOM6':'saddlebrown',
'NEMO3.6': 'green',
'Nemo3.4.1': 'lime',
'NEMO-HadGEM3-GO6.0':'darkturquoise',
'NEMO-LIM3.3.6':'darkkhaki',
'GISS Ocean':'teal',
'POP2': 'lightcoral',
'MPIOM1.6.3':'purple',
'MRI.COM4.4':'lightskyblue',
'MPAS-Ocean':'slategrey',
'EC-Earth NEMO3.6': 'blue',
}
# Create a new column 'Color' in the model_info_df DataFrame based on the ocean components
model_info_df['Color'] = model_info_df['Ocean Component'].map(ocean_component_colors)
# Rest of your code
input_folder = '/Volumes/Seagate Hub/cmip6_downloader/data/Salinity models/'
myfiles = os.listdir(input_folder)
for f in myfiles:
if '._so' in f:
myfiles.remove(f)
fig, axs = plt.subplots(nrows=1, ncols=2, figsize=(29, 12), sharex=True)
plt.subplots_adjust(right=0.7, hspace=0.3)
def plot_salinity_profile_OC(ax, file, color, start_year, end_year):
ds = xr.open_dataset(input_folder + file)
if 'CESM' in file:
ds['lev'] = ds['lev'] / 100
model_name = file.split('so.')[-1].split('.')[0]
ds_mean = ds.so.isel(time=np.arange((start_year-1850)*12, (end_year+1-1850)*12, 1)).mean(dim='time')
y = ds_mean.mean(dim=['longitude', 'latitude'])
ax.plot(y, ds['lev'], label=model_name, color=color, linewidth=2, linestyle='-')
# Rest of the plot_salinity_profile function remains the same
for file in myfiles:
model_name = file.split('so.')[-1].split('.')[0]
color = model_info_df.loc[model_info_df['Model Name'] == model_name, 'Color'].values[0]
plot_salinity_profile_OC(axs[0], file, color, start_year=1970, end_year=1980)
plot_salinity_profile_OC(axs[1], file, color, start_year=2000, end_year=2010)
axs[0].plot(AJX_avg_salinity, depth_values, label='1975 AIDJEX', color='red', linewidth = 5)
axs[1].plot(ITP_avg_salinity, depth_values, label='2006-2012 ITP', color='orange', linewidth=5)
for i, ax in enumerate(axs):
ax.set_ylabel('Depth (m)', fontsize=14)
ax.set_xlabel('Salinity (g/kg)', fontsize=14)
ax.set_ylim(0, 150)
ax.set_xlim(27, 35)
ax.tick_params(axis='y', labelsize=12)
ax.tick_params(axis='x', labelsize=12)
ax.invert_yaxis()
label = chr(99 + i) # Convert index to ASCII character ('c', 'd', etc.)
ax.text(-0.1, 1, label, fontsize=19, fontweight='bold', transform=ax.transAxes)
# ... rest of your code above
unique_ocean_components = model_info_df['Ocean Component'].dropna().unique()
unique_colors = [ocean_component_colors[oc] for oc in unique_ocean_components]
# Create custom legend handles
legend_handles = [plt.Line2D([0], [0], color=color, lw=2, label=ocean_component) for ocean_component, color in zip(unique_ocean_components, unique_colors)]
# Add custom legend handles for '1975 AIDJEX' and '2006-2012 ITP'
aidjex_handle = plt.Line2D([0], [0], color='red', lw=5, label='1975 AIDJEX')
itp_handle = plt.Line2D([0], [0], color='orange', lw=5, label='2006-2012 ITP')
legend_handles.extend([aidjex_handle, itp_handle])
# Create a shared legend for all subplots
fig.legend(legend_handles, [handle.get_label() for handle in legend_handles], fontsize=15, bbox_to_anchor=(0.68, -0.1, 0.15, 1), ncol=1, borderaxespad=0)
plt.savefig('ocean_component_plot.png', bbox_inches='tight')
# Read the Excel sheet into a pandas DataFrame
excel_file_path = '/Volumes/Seagate Hub/ocean_components.xlsx'
model_info_df = pd.read_excel(excel_file_path, engine='openpyxl', sheet_name='Sheet1', header=None, names=['Model Name', 'Ocean Component'])
input_folder = '/Volumes/Seagate Hub/cmip6_downloader/data/Salinity models/'
myfiles = os.listdir(input_folder)
for f in myfiles:
if '._so' in f:
myfiles.remove(f)
fig, axs = plt.subplots(nrows=1, ncols=2, figsize=(29, 12), sharex=True)
plt.subplots_adjust(right=0.7, hspace=0.3)
unique_ocean_components = model_info_df['Ocean Component'].dropna().unique()
color_palette = sns.color_palette('tab10', n_colors=len(unique_ocean_components))
ocean_component_colors = dict(zip(unique_ocean_components, color_palette))
# Create a new column 'Color' in the model_info_df DataFrame based on the ocean components
model_info_df['Color'] = model_info_df['Ocean Component'].map(ocean_component_colors)
# Rest of your code
def plot_salinity_profile_OC(ax, file, color, start_year, end_year):
ds = xr.open_dataset(input_folder + file)
if 'CESM' in file:
ds['lev'] = ds['lev'] / 100
model_name = file.split('so.')[-1].split('.')[0]
ds_mean = ds.so.isel(time=np.arange((start_year-1850)*12, (end_year+1-1850)*12, 1)).mean(dim='time')
y = ds_mean.mean(dim=['longitude', 'latitude'])
ax.plot(y, ds['lev'], label=model_name, color=color, linewidth=2, linestyle='-')
# Rest of the plot_salinity_profile function remains the same
for file in myfiles:
model_name = file.split('so.')[-1].split('.')[0]
color = model_info_df.loc[model_info_df['Model Name'] == model_name, 'Color'].values[0]
plot_salinity_profile_OC(axs[0], file, color, start_year=1970, end_year=1980)
plot_salinity_profile_OC(axs[1], file, color, start_year=2000, end_year=2010)
axs[0].plot(AJX_avg_salinity, depth_values, label='1975 AIDJEX', color='red', linewidth = 5)
axs[1].plot(ITP_avg_salinity, depth_values, label='2006-2012 ITP', color='orange', linewidth=5)
for i, ax in enumerate(axs):
ax.set_ylabel('Depth (m)', fontsize=14)
ax.set_xlabel('Salinity (g/kg)', fontsize=14)
ax.set_ylim(0, 150)
ax.set_xlim(27, 35)
ax.tick_params(axis='y', labelsize=12)
ax.tick_params(axis='x', labelsize=12)
ax.invert_yaxis()
label = chr(99 + i) # Convert index to ASCII character ('c', 'd', etc.)
ax.text(-0.1, 1, label, fontsize=19, fontweight='bold', transform=ax.transAxes)
# ... rest of your code above
unique_ocean_components = model_info_df['Ocean Component'].dropna().unique()
unique_colors = [ocean_component_colors[oc] for oc in unique_ocean_components]
# Create custom legend handles
legend_handles = [plt.Line2D([0], [0], color=ocean_component_colors[oc], lw=2, label=oc) for oc in unique_ocean_components]
aidjex_handle = plt.Line2D([0], [0], color='red', lw=5, label='1975 AIDJEX')
itp_handle = plt.Line2D([0], [0], color='orange', lw=5, label='2006-2012 ITP')
legend_handles.extend([aidjex_handle, itp_handle])
# Create a shared legend for all subplots
fig.legend(legend_handles, [handle.get_label() for handle in legend_handles], fontsize=15, bbox_to_anchor=(0.68, -0.1, 0.15, 1), ncol=1, borderaxespad=0)
plt.savefig('ocean_component_plot.png', bbox_inches='tight')
import os
import numpy as np
import pandas as pd
import xarray as xr
import seaborn as sns
import matplotlib.pyplot as plt
input_folder = '/Volumes/Seagate Hub/cmip6_downloader/data/Salinity models/'
myfiles = os.listdir(input_folder)
for f in myfiles:
if '._so' in f:
myfiles.remove(f)
colors = sns.color_palette("husl", len(myfiles))
linestyles = ['-', '--', '-.', ':'] * (len(myfiles) // 4 + 1)
fig, axs = plt.subplots(nrows=2, ncols=2, figsize=(25, 20), sharex=True)
plt.subplots_adjust(right=0.7, hspace=0.3)
def plot_salinity_profile(ax, file, color, linestyle, start_year, end_year):
ds = xr.open_dataset(input_folder + file)
if 'CESM' in file:
ds['lev'] = ds['lev'] / 100
model_name = file.split('so.')[-1].split('.')[0]
ds_mean = ds.so.isel(time=np.arange((start_year-1850)*12, (end_year+1-1850)*12, 1)).mean(dim='time')
y = ds_mean.mean(dim=['longitude', 'latitude'])
ax.plot(y, ds['lev'], label=model_name, color=color, linewidth=2, linestyle=linestyle)
for file, color, linestyle in zip(myfiles, colors, linestyles):
plot_salinity_profile(axs[0, 0], file, color, linestyle, 1970, 1980)
plot_salinity_profile(axs[0, 1], file, color, linestyle, 2000, 2010)
axs[0, 0].plot(AJX_avg_salinity, depth_values, label='1975 AIDJEX', color='red', linewidth=4)
axs[0, 1].plot(ITP_avg_salinity, depth_values, label='2006-2012 ITP', color='orange', linewidth=4)
for i in range(2):
for j in range(2):
ax = axs[i, j]
ax.set_ylabel('Depth (m)', fontsize=14)
ax.set_xlabel('Salinity (g/kg)', fontsize=14)
ax.set_ylim(0, 150)
ax.set_xlim(27, 35)
ax.tick_params(axis='y', labelsize=12)
ax.tick_params(axis='x', labelsize=12)
ax.invert_yaxis()
# Add subplot labels
label = chr(97 + i * 2 + j) # Convert index to ASCII character ('a', 'b', etc.)
ax.text(-0.1, 1, label, fontsize=19, fontweight='bold', transform=ax.transAxes)
axs[0, 0].set_title('1970-1980 Salinity Profiles', fontsize=20)
axs[0, 1].set_title('2000-2010 Salinity Profiles', fontsize=20)
def plot_salinity_profile_ocean_component(ax, file, color, start_year, end_year):
ds = xr.open_dataset(input_folder + file)
if 'CESM' in file:
ds['lev'] = ds['lev'] / 100
model_name = file.split('so.')[-1].split('.')[0]
ds_mean = ds.so.isel(time=np.arange((start_year-1850)*12, (end_year+1-1850)*12, 1)).mean(dim='time')
y = ds_mean.mean(dim=['longitude', 'latitude'])
ax.plot(y, ds['lev'], label=model_name, color=color, linewidth=2, linestyle='-')
for file in myfiles:
model_name = file.split('so.')[-1].split('.')[0]
color = model_info_df.loc[model_info_df['Model Name'] == model_name, 'Color'].values[0]
plot_salinity_profile_ocean_component(axs[1, 0], file, color, start_year=1970, end_year=1980)
plot_salinity_profile_ocean_component(axs[1, 1], file, color, start_year=2000, end_year=2010)
axs[1, 0].plot(AJX_avg_salinity, depth_values, label='1975 AIDJEX', color='black', linewidth=4)
axs[1, 1].plot(ITP_avg_salinity, depth_values, label='2006-2012 ITP', color='black', linewidth=4)
# Get the handles and labels from the first subplot
handles, labels = axs[0, 0].get_legend_handles_labels()
# Add the handle and label for the ITP line from the second subplot
handles_itp, labels_itp = axs
import os
import xarray as xr
import numpy as np
import matplotlib.pyplot as plt
input_folder = '/Volumes/Seagate Hub/cmip6_downloader/data/Salinity models/'
myfiles = [f for f in os.listdir(input_folder) if not f.startswith('._') and 'so' in f]
fig, axs = plt.subplots(8, 4, layout='tight', figsize=(5 * 5, 5 * 8), sharex=True)
axs[-1, -1].remove()
axs[-1, -2].remove()
axs = axs.ravel()
def int_to_roman(integer):
val = [
1000, 900, 500, 400,
100, 90, 50, 40,
10, 9, 5, 4,
1
]
syb = [
"M", "CM", "D", "CD",
"C", "XC", "L", "XL",
"X", "IX", "V", "IV",
"I"
]
roman_num = ''
i = 0
while integer > 0:
for _ in range(integer // val[i]):
roman_num += syb[i]
integer -= val[i]
i += 1
return roman_num
# Extract salinity at depth 0 for each month from May to December
SA_depth0 = SA[4:12, 0] # May is at index 4 (0-based), depth 0 is at index 0
SI_depth0 = SI[4:12, 0]
# Create an array for the months (May to December)
months_obs = ['May', 'Jun', 'Jul', 'Aug', 'Sep', 'Oct', 'Nov', 'Dec']
def monthly_mean(file, i):
print(file)
ds = xr.open_dataset(input_folder + file)
if 'CESM' in file:
ds['lev'] = ds['lev'] / 100
model_name = file.split('so.')[-1].split('.')[0]
ds_month = ds.resample(time="M").mean()
ds_month_mean = ds_month.mean(dim=['longitude', 'latitude'])
ds_month_mean_depthzero = ds_month_mean.isel(lev=0)
months_to_select = range(5, 13) # May to December
decades = [(1970, 1980), (1980, 1990), (1990, 2000), (2000, 2010)]
for j, (start_year, end_year) in enumerate(decades):
month_indices = [(int(year) - 1850) * 12 + month - 1 for year in range(start_year, end_year) for month in months_to_select]
monthly_decadal_means = ds_month_mean_depthzero['so'].isel(time=month_indices).groupby('time.month').mean(dim='time')
monthly_decadal_std = ds_month_mean_depthzero['so'].isel(time=month_indices).groupby('time.month').std(dim='time')
y_upper = monthly_decadal_means + monthly_decadal_std
y_lower = monthly_decadal_means - monthly_decadal_std
axs[i].plot(months_to_select, monthly_decadal_means, label=f"{start_year}-{end_year}", linewidth = 4, color=plt.cm.viridis(j/len(decades)))
axs[i].fill_between(months_to_select, y_upper, y_lower, color=plt.cm.viridis(j/len(decades)), alpha=0.2)
axs[i].set_xlabel('Month', fontsize= 20)
axs[i].set_ylabel('Salinity (g/kg)', fontsize= 20)
axs[i].set_title(model_name, fontsize = 25)
axs[i].set_ylim(25,33)
axs[i].plot(months_to_select, SA_depth0, label='AIDJEX', color = 'red', linewidth = 4)
axs[i].plot(months_to_select, SI_depth0, label='ITP', color = 'orange', linewidth =4)
if i == len(axs) - 3:
axs[i].legend(loc='lower left')
axs[i].tick_params(axis='y', labelsize=18)
axs[i].tick_params(axis='x', labelsize=18)
axs[i].xaxis.set_tick_params(which='both', labelbottom=True)
if i < 26:
subplot_label = string.ascii_lowercase[i]
else:
subplot_label = int_to_roman((i - 25))
axs.reshape(-1)[i].text(0.02, 1.1, subplot_label, transform=axs.reshape(-1)[i].transAxes, fontsize=16, fontweight='bold', va='top', ha='left')
i = 0
for file in myfiles:
monthly_mean(file, i)
i += 1
plt.savefig('figure2.png')
so.BCC-ESM1.CadBasin.nc so.CNRM-CM6-1-HR.CadBasin.nc so.CAMS-CSM1-0.CadBasin.nc so.GISS-E2-1-G.CadBasin.nc so.CanESM5.CadBasin.nc so.CNRM-CM6-1.CadBasin.nc so.CESM2-WACCM-FV2.CadBasin.nc so.TaiESM1.CadBasin.nc so.CESM2.CadBasin.nc so.CNRM-ESM2-1.CadBasin.nc so.CESM2-WACCM.CadBasin.nc so.MPI-ESM1-2-HR.CadBasin.nc so.BCC-CSM2-MR.CadBasin.nc so.CESM2-FV2.CadBasin.nc so.EC-Earth3-CC.CadBasin.nc so.EC-Earth3-Veg-LR.CadBasin.nc so.EC-Earth3.CadBasin.nc so.E3SM-1-0.CadBasin.nc so.HadGEM3-GC31-LL.CadBasin.nc so.CMCC-CM2-SR5.CadBasin.nc so.CMCC-ESM2.CadBasin.nc so.EC-Earth3-AerChem.CadBasin.nc so.EC-Earth3-Veg.CadBasin.nc so.FIO-ESM-2-0.CadBasin.nc so.GFDL-CM4.CadBasin.nc so.GFDL-ESM4.CadBasin.nc so.GISS-E2-2-G.CadBasin.nc so.MPI-ESM-1-2-HAM.CadBasin.nc so.MPI-ESM1-2-LR.CadBasin.nc so.MRI-ESM2-0.CadBasin.nc